DEVELOPMENT  OF  A  DEMAND  FORECASTING  MODEL 

FOR 

FLORAL  DATA 


by  i"4^ 

WILLIAM  FRANCIS  IRVINE 
B.  S.',  University  of  Notre  Dame,  1966 


A  MASTER'S  THESIS 

submitted  in  partial  fulfillment  of  the 

^    requirements  for  the  degree 

MASTER  OF  SCIENCE 

Department  of  Industrial  Engineering 


KANSAS  STATE  UNIVERSITY 
Manhattan,  Kansas 


1968 


Approved  by: 


J-  C-  J{j/t^<ii^ 


Major  Professor 


LD  ■ 

Ti 

I  US' 

c  .o- 


TABLE  OF  CONTENTS 


Page 
1 


INTRODUCTION  

CHAPTER  I.   Representation  of  the  Time  Series 5 

CHAPTER  II.   Data  Analysis 16 

CHAPTER  III.   Forecasting  the  Time  Series 45 

CHAPTER  IV.   Choosing  the  Smoothing  Constant  and  the  Initial 

Coefficient  Vector   65 

CHAPTER  V.   Gathering  the  Data 80 

CHAPTER  VI.   Discussion  of  Results 86 

CHAPTER  VII.   Conclusions 106 

BIBLIOGRAPHY 109 

ACKHOWLEDGEiffiNTS 110 


APPENDIX  A 
APPENDIX  B 
APPENDIX  C 
APPENDIX  D 
APPENDIX  E 
APPENDIX  F 


Derivation  of  Statistical  Regression  Equation  Ill 

Derivation  of  Exponential  Smoothing  Equation 114 

Polynomial  Regression 118 

The  Z-Transform 120 

Tables  and  Figures 122 

Computer  Programs 184 


INTRODUCTION 

A  common  problem  in  business  and  industry  today  is  the  accurate 
estimation  of  the  level  of  inventory  necessary  to  meet  some  future 
demand.   In  some  situations  this  task  may  be  quite  critical.   This  can 
best  be  explained  using  an  example.   Since  this  research  concerns  the 
floral  industry,  let  us  use  as  an  example  the  demand  for  a  floral  product. 

Consider  a  plant  which  has  a  growing  period  of  two  months  from  the 
time  it  is  first  planted  until  it  is  ready  to  be  sold.   In  order  to 
know  how  many  plants  to  start  now,  we  must  have  some  idea  of  how  many 
will  be  needed  for  sale  two  months  from  now.   If  we  underestimate  this 
future  demand,  we  will  not  have  sufficient  plants  in  stock  that  are 
ready  for  sale.   Conversely,  if  we  overestimate  future  demand  we  will 
have  too  many  plants  and  may  be  forced  to  sell  them  at  a  loss  or  even 
throw  them  away.   Thus  we  see  that  it  is  desirable  to  be  able  to  fore- 
cast this  future  demand  with  the  confidence  that  this  forecast  will  not 
be  too  far  from  the  actual  demand  two  months  after  planting.   The  impor- 
tance of  forecasting  can  be  seen  almost  anywhere  an  inventory  is  main- 
tained. 

The  logical  basis  for  a  forecast  of  demand  for  any  item  is  the 
demand  history  of  that  item.   Normally,  such  past  demand  data  is  in  the 
form  of  a  time  series,  a  sequence  of  observations  made  at  specific  and 
equal  time  intervals.   There  are  several  ways  in  which  future  demands 
may  be  interpretted  from  this  time  series.  Among  these  methods,  subject- 
ive reasoning  is  probably  the  simplest.  This  is  merely  the  opinion  of 


an  experienced  manager  of  what  demand  will  be.  Other  methods  include 
statistical  regression  and  graphical  curve  fitting.   All  of  these 
methods  nay  be  classified  as  statLc;  that  is,  they  deal  in  some  way  with 
fitting  a  curve  to  the  data  as  it  now  exists.  Perhaps  a  more  realistic 
approach  is  to  give  the  recent  data  more  weight  than  older  data  on  the 
assumption  that  this  newer  data  more  accurately  foretells  the  behavior  of 
demand  in  the  future.  Methods  employing  such  weighting  factors  may  be 
classified  as  dynamic.  The  moving  average  is  one  approach  to  weighting. 
This  method  takes  an  average  of  the  last  N  data  values  as  a  forecast  of 
demand  for  the  next  sampling  interval,  where  N  is  some  pre-determlned 
number  of  sampling  intervals.  Notice  that  if  N  is  large  we  consider 
much  past  data,  and  if  N  is  small  we  put  the  emphasis  on  only  the  most 
recent  data.  Thus  we  may  control  the  rate  of  response  of  this  forecast- 
ing system  by  varying  N.   If  we  wish  to  include  all  past  data,  we  may  use 
the  general  exponential  smoothing  technique.  This  is  the  method  used  for 
this  research.  General  exponential  smoothing  gives  data  taken  k  periods 
in  the  past  the  weightaSk^  where  a  and B  are  constants  and  their  sum  equals 
one.  Now  we  may  control  the  system  response  rate  to  changes  in  the  data 
by  varying  the  value  of  B,  known  as  the  discount  factor. 

General  exponential  smoothing  essentially  fits  a  least  squares  curve 
to  the  weighted  time  series.   By  least  squares  we  mean  that  the  curve 
is  determined  In  a  manner  that  minimizes  the  sum  of  weighted  square 
residuals,  or  fitting  errors.  This  method  provides  added  adaptability 
to  changes  in  the  data  by  re-evaluating  the  model  coefficients  at  each 
sampling  interval.  Let  us  illustrate  the  meaning  of  this  statement.  We 


may  describe  the  time  series  by: 

x(t)  =  5(t)  +  e(t) 

where     x(c)  =  the  past  data  time  series  at  time  t 

£;(t)  =  the  systematic  component  of  x(t)  at  time  t 

e(t)  =  some  random  variation  of  x(t)  from  5(t)  at  time  t. 

We  are  assuming,  of  course,  that  the  time  series  contains  some  mathemat- 
ically describable  process  ^(t).   We  wish  to  formulate  a  model  to  describe 
C(t)  of  the  form 


x(t)  -  E  a  (t)  f  (t) 

i=l 


where     x(t)  =  C(t),  if  the  model  above  is  the  correct  one 
f  (t)  =  a  deterministic  function  of  time 

=  the  i''^  fitting  function 
a^(t)  =  the  least  squares  coefficient  of  f.(t). 

Now,  by  a  re-evaluation  of  model  fit  we  mean  that  general  exponential 
smoothing  re-computes  the  least  squares  fitting  function  coefficients 
of  the  model  at  each  sampling  interval.   In  this  updating  procedure  the 
error  made  on  the  previous  forecast  is  taken  into  account.   Thus,  we 
may  note  several  advantages  of  general  exponential  smoothing  over  the 
methods  discussed  earlier.   First,  it  accounts  for  all  past  data. 
Secondly,  exponential  smoothing  gives  the  investigator  more  immediate 
control  over  the  forecasting  system  response  rate.  As  is  noted  in  the 
text  to  follow,  this  method  also  requires  less  data  storage  once  it  is 


started,  since  all  past  data  contained  in  one  word  of  information.   A 
fourth  advantage  of  general  exponential  smoothing  is  the  re-evaluation 
of  the  least  squares  coefficients  in  the  light  of  each  additional  amount 
of  data. 

The  methods  used  here  are  taken  primarily  from  Smoothing,  Forecasting 
and  Prediction  of  Discrete  Time  Series  by  Robert  Goodell  Brown  (2) .   All 
computations  are  performed  on  the  IBM  System  360/50  digital  computer  at 
Kansas  State  University.   The  data  used  consists  of  four  discrete  time 
series  of  demand  for  potted  plants: 

(1)  Wholesale  demand  for  chrysanthemums  in  monthly  dollar  receipts. 

(2)  Retail  demand  for  chrysanthemums  in  monthly  dollar  receipts. 

(3)  Wholesale  demand  for  lilies  in  daily  plant  sales. 

(4)  Wholesale  demand  for  poinsettias  in  daily  plant  sales. 

All  time  series  were  taken  from  the  sales  records  of  Manhattan  Retail 
and  Wholesale  Floral  Companies  in  Manhattan,  Kansas.   Chrysanthemums 
are  sold  throughout  the  entire  year,  while  lilies  and  poinsettias  are 
seasonal  items,  sold  only  during  the  Easter  and  Christmas  seasons, 
respectively. 

This  particular  data  was  used  due  to  the  interest  of  the  author. 
Having  been  employed  by  Manhattan  Floral  Company  for  several  years,  the 
author  feels  that  he  has  a  greater  understanding  of  the  problems  involved 
in  this  area  than  in  any  other.   In  addition,  data  from  this  source 
could  be  more  easily  obtained  and  verified. 


1.   REPRESENTATION  OF  THE  TIME  SERIES 

1.1  Form  of  the  Time  Series  Model 

A  time  series  is  a  set  of  observations  taken  at  specific  points  in 
time,  usually  at  equal  intervals.   The  analyses  that  follow  assume  that 
the  time  series  involved  may  be  separated  into  two  components,  systematic 
and  random.   The  systematic  component  can  be  represented  as  a  determinis- 
tic function  of  time.   Thus,  the  time  series  may  be  represented  mathemat- 
ically as: 


^  =  ^t  +  ^t 


where 

X  =  the  actual  time  series  observation  at  time  t. 

E,^  =  the  systematic  component  of  the  observation  at  time  t. 

Ej^  =  the  random  component  of  the  observations  at  time  t. 

The  distribution  of  the  random  component  of  the  time  series,  Ej.,  possesses 
the  following  properties: 


E[ej]  -  Expected  value  of  e^  -  0 

E[Ej  •  c^_j]  =  0  for  ji^O 

-  0^  for  j-0 

where  a  is  the  variance  of  the  random  noise  distribution.   Successive 
samples  of  the  random  component  are  assumed  to  have  no  serial  correlation. 
Having  defined  e^  as  random  with  a  mean  value  of  zero,  we  make  no 


attempt  to  describe  it.   So,  in  order  to  represent  the  time  series  x  , 
except  for  the  random  variation,  e  ,  we  must  be  able  to  describe  the 
systematic  component,  Cj.   Once  we  are  able  to  describe  C^,  we  may  pre- 
dict future  values  of  Xj^,  since  ?(.  is  a  function  of  time  only  and  time 
Is  always  known  for  the  future. 

To  better  visualize  the  time  series,  let  us  look  at  the  plot  shown 
below. 


x(t) 


• 

..---' 

,^y 

/^-^     past   time 

future  time 

Fig.  1:1:1  Plot  of  hypothetical  time  series 


This  is  the  plot  of  a  time  aeries  x  on  the  vertical  axis,  versus 
time  on  the  horizontal  axis.   Xj.  may  be  any  time  series  for  which  we 
wish  to  predict  future  values. 

The  following  are  a  few  of  the  methods  that  may  be  employed  to 
attack  this  forecasting  problem: 

1.  Subjective  estimate 

2.  Graphical  curve  fitting 

3.  Statistical  regression  analysis 
A.   Simulation 

5 .  Moving  averages 

6.  Exponential  smoothing. 


For  reasons  that  will  be  explained  in  sections  1.2  and  1.3,  we  may 
classify  the  first  three  methods  above  as  static  methods,  and  the  last 
three  as  dynamic  methods. 

1.2  Static  Methods  of  Time  Series  Representation 

The  first  of  the  static  methods,  subjective  estimate,  implies  an 
opinion,  guided  by  reason,  experience,  and  judgment,  of  what  x^.  will  be 
for  a  particular  point  in  time,  t.  This  prediction  would  normally  be 
made  by  a  person  or  persons  familiar  with  the  time  series  involved.   It 
is  not  based  on  any  particular  analytical  foundation. 

Graphical  curve  fitting  may  be  accomplished  by  merely  observing  the 
time  series  plot,  like  the  one  in  Figure  1:1:1,  and  drawing  a  representa- 
tive curve  through  the  data  points  shown  in  the  "past  time"  section  of 
the  plot.  These  points  are  actual  data  gathered  from  past  performance 
of  X.  . 

By  extending  this  curve  into  "future  time"  a  prediction  of  the 
future  values  of  x^^  may  be  made.   This  method,  though  not  totally  subjec- 
tive, may  fail  to  adequately  represent  any  cyclic  fluctuation  inherent 
in  the  time  series. 

Probably  the  major  analytical  method  of  static  time  series  analysis 
is  statistical  regression.   This  method  attempts  to  formulate  a  model 
that  will  describe  the  systematic  component  of  the  time  series,  5^..   The 
model  may  contain  any  independent  variables  that  are  known  exactly  for 
all  future  time.   It  is  up  to  the  investigator  to  choose  the  proper  terms 
and  their  form  within  the  model.  Statistical  regression  techniques  are 


used  to  determine  the  proper  coefficients  of  the  terms  in  the  model. 
These  techniques  will  give  coefficients  that  effect  a  minimization  of 
the  sum  of  squared  residuals,   S,  for  the  time  series  values  given 
(Appendix  A) . 

S  -  E   [x  -  X  ] 

t-0 


N  =  the  number  of  data  values  given 
Xf  "  the  actual  time  series  value  for  time  t 
X  =  the  time  series  value  for  time  t  as  determined  from 
the  model. 

Although  this  method  minimizes  the  sum  of  the  squared  fitting  errors  for 
the  data  given,  it  makes  no  provision  for  a  change  in  the  general  form 
of  the  time  series  through  time.   As  is  true  with  all  static  methods  of 
time  series  representation,  statistical  regression  makes  no  attempt  to 
reevaluate  the  model  coefficients  at  each  point  in  time  to  reflect  the 
forecast  error  made  in  the  previous  prediction. 

1.3  Dynamic  Methods  of  Time  Series  Representation 

Simulation  techniques  attempt  to  formulate  a  process  model  that  will 
behave  like  the  actual  process  over  some  specified  period  of  time.   For 
example,  we  attempt  to  establish  probability  distributions  that  are 
representative  of  the  variables  contributing  to  the  process,  such  as 
demand,  production,  etc.   These  distributions  will  generally  be  based  on 


past  behaviour  of  the  actual  process,  if  such  data  is  available.   By 
using  a  model  to  simulate  the  true  process  over  some  period  of  time,  we 
may  gain  an  insight  as  to  what  values  the  variables  will  exhibit  in  the 
future. 

Normally  we  would  expect  the  more  recent  values  in  the  time  series 
to  be  more  representative  of  future  values  than  the  older  data.   The  last 
two  dynamic  methods  to  be  discussed  account  for  this  probability  by  paying 
much  attention  to  recent  data  and  little  or  no  attention  to  older  data. 
This  important  feature  of  dynamic  methods  of  time  series  representation 
allows  the  estimates,  x  to  reflect  a  change  in  the  form  of  the  time 
series  as  time  progresses. 

The  moving  average  is  a  simplified  form  of  dynamic  time  series  repre- 
sentation. For  purposes  of  illustration,  consider  a  time  series  for  which 
the  systematic  component  may  be  represented  as  a  constant  process. 


where 

x^  =  the  observed  time  series  value  for  time  t 

a^   =  5(t)  =  the  value  of  a  constant  process  at  time  t 

Ej  =  the  random  variation  from  the  true  process  at  time  t. 

It  will  become  clear  in  the  explanation  that  follows  that  a^.  may  actually 
change  slowly  with  time.   Since  this  is  a  dynamic  representation,  this 
change  in  the  true  process  will  be  accounted  for  in  the  analysis. 

Since  we  have  defined  the  systematic  component  of  the  time  series, 
5j,  as  a  deterministic  function  of  time,  we  may  construct  a  function,  or 
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model,  that  will  describe  the  time  series,  except  for  a  random  variation, 
5^.   In  this  case  a  is  our  model  of  the  true  process,  with  time  as  the 
only  independent  variable. 

As  an  estimate  of  ap^.^^,  we  might  take  an  average  of  all  the  observa- 
tions, X  ,  up  to  the  present  time,  T.   So, 


t=0 


where 


a    =  the  estimated  value  of  the  true  process,  a  ,  at 
T+1  "^ 

time  t  =  T  +  1. 


A  better  estimate  of  a„^,  might  be  an  average  over  only  the  most  recent 
N  observations.   This  is  what  is  known  as  a  moving  average  M^.  Now, 


t+1 


■'t-1  ^  -^t  ~  "^t-N 


M  =  the  moving  average  of  x  over  the  most  recent 
N  observations,  or  the  moving  average  of 
period  N. 


M^  minimizes  the  sum  of  squares  of  differences  between  the  most  recent 
N  observations  and  our  estimates  a  of  the  process  model.  Note  that 
when  N  is  large  we  are  giving  more  weight  to  older  data  than  when  N  is 
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small. 

The  moving  average,  though  it  is  a  dynamic  method  of  forecasting, 
still  has  two  important  limitations.   First,  it  does  not  consider  all 
past  data,  and  gives  equal  weight  to  the  observations  that  it  does  con- 
sider.  Second,  it  is  difficult  to  change  the  response  rate  of  a  system 
based  on  moving  averages.   Though  the  system  response  is  based  on  the 
value  of  N,  it  would  be  difficult  to  alter  a  record  system  based  on  a 
particular  N  value.   Also,  if  there  are  significant  changes  in  the  true 
process  through  time,  the  moving  average  will  give  highly  unsatisfactory 
forecasts. 

Again  for  purposes  of  illustration,  suppose  that  we  lost  all  our 
data,  except  the  value  of  the  previous  moving  average  M.j._-,  •  Then,  by 
this  same  scheme,  our  best  estimate  for  x^^^.   would  be, 

N 


ix.j,  +  (I  -  |)  Mj,.^ 


where 


H-  =  an  estimate  of  M„. 


T 
This  process  is  known  as  "smoothing."  Now  let, 


1 

"  =  N 

S^(x)  -  Mj.. 
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We  now  have 


where 


S^(x)  =  axj.  +  3S^_^  (x) 


S|.(x)  =  the  smoothed  value  of  x  at  time  t. 

S  , ,  N  =  the  smoothed  value  of  x  at  time  t-1. 
t-l(x) 

a  =  the  smoothing  constant 
3  =  the  discount  factor 

a  and  B  are  not  exactly  the  same  as  —  and  (1  -  i) ,  tut  are  similar  to 
them.   We  shall  see  their  significance  later  in  this  chapter.   Note  here 
that. 


a,  g  ;  1 
a  +  B  =  1. 

This  operation,  carried  out  over  a  sequence  of  observations,  is  called 
exponential  smoothing.   To  see  this  effect  going  back  in  time  we  merely 
expand  our  original  smoothing  equation. 


S^(x)  =  axj.  +  6S^_^(x) 

=  axj.  +  3[ax^_^  +  6S^_2(x)] 


2 
ax  +  ciBx  j^  +  a6  x 


+  B  X  . 
o 
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The  exponent  of  B  increases  geometrically  as  we  go  back  in  time.  So,  the 
weight  given  the  initial  x  value,  x  ,  after  k  observations  is  8  .  Hence, 
we  call  this  exponential  smoothing. 

This  procedure  offers  several  advantages  over  the  methods  discussed 
earlier.   First,  all  the  past  data  is  taken  into  consideration.   The 
older  the  data  is,  the  less  weight  it  receives  when  determining  Sj.(x). 
Also,  we  may  change  the  response  rate  of  the  system  by  merely  altering 
a,  the  smoothing  constant,  or  B,  the  discount  factor.   So,  for  B  small, 
we  are  giving  more  weight  to  past  data  than  if  B  were  larger.   If  6  is 
large,  then,  the  system  will  adjust  more  quickly  to  a  wide  fluctuation 
in  the  observations.   The  third  advantage  of  smoothing  is  that  all  past 
data  is  presented  by  only  one  value,  S(.(x),  instead  of  N  observations 
as  in  the  case  of  the  moving  average.   Lastly,  exponential  smoothing 
minimizes  the  weighted  sum  of  squared  residuals  (Appendix  B) . 

So  far  in  our  discussion,  we  have  been  dealing  with  the  simplest 
of  processes,  the  constant  process.   Let  us  now  consider  a  general  model 
of  the  true  process. 

x(t)  =  a^(t)  f^(t)  +  a^Ct)  f2(t)  +  .  .  .  +  a^(t)f^(t)  +  e(t) 

n 
=  I     a,(t)  f  (t)  +  c(t) 
1-1 


where 


th 
,  (t)  =  the  coefficient  of  the  i^^h  fitting  function  at  time  t. 


f  (t)  =  the  i"-"  fitting  function  at  time  t 


See  section  2.1  for  a  discussion  of  allowable  fitting  functions.   In  this 
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expression  the  systematic  component  of  the  time  series  Is  represented  by 


I      a.(t)  f  (t), 

1=1 


which  then  represents  the  model  of  our  system.   The  fitting  functions  are 
deterministic  functions  of  time.   In  the  case  of  the  constant  process, 
we  had  only  one  fitting  function  whose  value  was  the  constant,  1.   The 
coefficients  a.(t)  are  treated  as  constants  over  any  local  period  of  time. 
Nevertheless,  the  underlying  process  5 (t)  can  be  going  through  a  slow 
random  walk  with  respect  to  one  or  more  of  the  coefficients.   As  will  be 
shown  later,  the  exponential  smoothing  procedure  accounts  for  this. 

For  future  reference,  let  us  write  the  general  time  series  represen- 
tation in  matrix  notation.  ' 


x(t)  =  I     a  (t)  f  (t)  +  e(t) 
1=1 


a'(t)  f(t)  +  (t) 


where 


a(t)  = 


a^(t) 
a^Ct) 


a_(t) 


,  and   f(t)  = 


f2(t) 


f  (O 


Hence 


x(t)  -  a'(t)  f(t). 
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It  should  be  apparent  at  this  point  that  the  first  requirement  in 
the  formulation  of  a  forecasting  procedure  for  a  time  series  Is  the 
proper  choice  of  the  fitting  functions  that  make  up  the  model.   There 
are  several  guidelines  we  may  follow  in  making  this  choice.   These  guide- 
lines are  discussed  in  the  next  chapter. 

The  exponential  smoothing  procedure  will  be  dealt  with  In  much 
greater  detail  In  Chapter  III. 
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2.   DATA  ANALYSIS 

2.1  Choice  of  Fitting  Functions  for  the  Model 

In  the  first  chapter  we  described  the  systematic  component  of  the 
time  series.   C   is  composed  of  fitting  functions,  f.(t),  and  their 
corresponding  coefficients,  a  (t) .   In  matrix  notation  this  is 

e^  =  i(t)  =  a'(t)  f{t) 

In  this  chapter  we  shall  discuss  some  methods  for  determining  the  proper 
fitting  functions  to  be  used  in  the  model,  x  . 

The  time  series  data  with  which  we  are  dealing  dictates  the  form 
of  the  model  fitting  functions.   The  problem  arises  in  determining  from 
the  data  those  fitting  functions  that  will  best  describe  the  process 
inherent  within  the  time  series.   A  discussion  of  the  possible  fitting 
function  forms  will  serve  to  clarify  the  alternatives  available.   Keep 
in  mind  that,  by  definition,  the  fitting  functions  must  be  known  exactly 
for  future  time. 

Brown  (2,58)  classifies  time  series  models  into  four  groups: 

1.  Algebraic  models 

2.  Trancendental  models 

3.  Composite  models 

4.  Regression  models. 

The  first  of  these  classifications,  algebraic  models,  takes  in  all 
polynomial  functions  of  time,  including  the  constant  model  with  which 
we  dealt  in  Chapter  1.   Our  forecast  in  the  case  of  a  constant  model, 
or  polynomial  model  of  degree  zero,  is 
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x(t:)  =  a(t)  f(t) 

=  a(t)  t° 

•=  a(t)  (1) 
x(t)  =  a(t) 

where  a(t)  =  an  estimate  from  the  smoothing  process  of  the  fitting  func- 
tion coefficient.   As  we  shall  show  in  Chapter  3,  the  coefficients  are 
re-estimated  by  the  exponential  smoothing  technique  each  time  additional 
data  has  become  available.   Thus  we  allow  for  gradual  changes  in  the 
process  described  by  our  model.   This  model  might  be  used  to  represent  a 
time  series  that  changes  gradually,  if  at  all,  from  a  constant  value.   A 
pictorial  example  of  this  type  of  process  is  shown  in  Figure  2:1:1.   The 
plot  is  not  shown  as  a  smooth  line  because  of  the  random  variation  contri- 
buted by  £(t),  the  "noise"  element  In  x(t) . 


x(t) 


Fig,  2:1:1  Constant  process  time  series 

Going  one  step  further,  consider  the  case  where  the  time  series 
exhibits  a  nearly  linear  form,  as  in  Figure  2:1:2.   A  first  degree  poly- 
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noraial  model 


x(t) 


Fig.  2:1:2  Linear  process  time  series 
would  probably  describe  this  process  adequately.   So, 


f^(t)  =  1 


and 


f^U)  =  t 


x(t)  =  I     a  At)   i  At) 
1-1 


aj^(t)  (1)  +  a^Ct)! 


where 


a(t) 


a(t)    f(t) 

"a^(t)~ 

a2(t) 

_                          _ 

and  f(t)  =  Ip 


Again,  the  plot  shows,  as  always,  the  existence  of  some  random  "noise" 
in  the  data. 

We  may  similarly  formulate  a  polynomial  model  of  any  desired  degree. 
Let  us  describe  an  n*"  degree  polynomial  model  as  one  containing  a  constant 
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term  and  a  polynomial  of  degree  n-1.   Then,  an  n   degree  polynomial 
model  would  be, 


x(t)  =  j:  a.  (t)  £  (t) 
i=l  ^     ^ 


a^(t)  (1)  +  n^Mt  +  a^Wt'- 
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a(l:)t 


n-1 


(n-1) 1 


where 


a(t)  fCc) 


a^(t) 


a(t) 


33(1) 


a_(t) 


and  f(t)  = 


1 
t/11 
t2/21 


t"-V(n-i): 


Note  that  the  (n-1)   degree  polynomial  model  requires  the  estimation  of 
n  fitting  function  coefficients. 

Brown  (2,  53)  notes  that  the  number  of  observations  in  the  span 
covered  by  the  model  should  be  several  times  the  number  of  degrees  of 
freedom,  n.   Logically,  forecasts  based  on  only  a  few  observations  would 
be  very  unreliable. 

The  question  now  arises  as  to  what  degree  of  polynomial  is  appropriate 
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for  a  particular  time  series.   If  the  observations  may  be  exactly  repre- 
sented by  a  polynomial  of  degree  n,  i.e.  no  "noise",  that  degree  may  be 
determined  by  successive  differences.   In  this  case, 


=  t  =  ° 


We  denote  a  first  difference  as 


and  a  second  difference  as 

=  ^(t)  -2^(t-l)+S^'-^> 
Thus,  a  difference  of  order  k  is  simply 

\  =  ^t  «  -  \-l^ 

If  the  appropriate  polynomial  is  of  degree  n,  the  difference  A    will 
be  identically  zero.   In  the  more  practical  case  in  which  e   is  not 
zero,  the  (n  +  1)  st  differences  of  the  data  A   x  will  average  zero 
over  time,  if  the  process  5  is  a  polynomial  of  degree  n.   The  variance 
of  the  differences  will  grow  as  the  sum  of  squares  of  the  binomial 
coefficients  in  the  expansion  of  A   x.  So,  if  e  has  a  variance  a  ,  the 
values  A   x  will  have  variance 

a^[2(n+l)]l/[(n+l)l]^ 
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2 
This  variance  is  larger  than  a  for  n+1  >  5  (Brown,  2,  63). 

In  many  cases  the  proper  polynomial  degree  may  be  determined  by 
merely  inspecting  a  plot  of  the  data  versus  time.   Statistical  least 
squares  regression  analysis  is  another  very  useful  method  of  attaining 
this  end.   This  last  approach  will  be  mentioned  in  our  discussion  of 
trend  analysis. 

Brown's  second  classification  of  time  series  models,  transcendental 
models,  is  characterized  by  two  sub-classifications: 

a.  Exponential  models 

b.  Trigonometric  models 

These  two  will  be  dealt  with  separately. 

An  exponential  function  may  be  used  to  describe  a  process  in  which 
the  rate  of  growth  is  proportional  to  the  state  of  growth  in  time,  like 
compound  interest.   Here,  the  change  in  value  from  5/  %  to  C(t+1)  may  be 
expressed  as  a  constant  percentage  or  C (t) ,  for  any  t.   The  form  of  5 
for  a  process  which  is  a  simple  exponential  function  is 


Note  that 


5(t-l)  =  ka''"^ 


at)/e(t-l)  =  kaVka 


t-1 


where 


a  constant  in  the  short  run. 
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This  simple  model,  when  plotted  on  semilogarithmic  paper,  will  define 
a  nearly  straight  line  with  slope  log  a,  since 


log  e.^)  =  log  k.  +  t  log 


The  general  form  of  the  exponential  model  is 


where 


5(t)  =  x(t)  =  k^(^)a'  +  k2(J)a'^"^b  +  Y^i^)  a!'~\'^ 


.  j_  1  /  t  ,  t-n+L  n-1 

+  .  .  .  +  k  (  ,)a    T) 
n  n-1 


V    (t-k)lkl 


is  the  binomial  coefficient  of  the  (k+l)st  term  in  the  expansion  of 
(a+b)^ 

To  describe  processes  which  are  periodic,  it  is  most  convenient  to 
use  the  second  of  the  trancendental  model  types,  trigonometric  models. 
By  periodic  we  mean  a  process  whose  time  plot  shows  repetitive  patterns. 
The  length  of  time  between  similar  sequence  patterns  is  called  the 
period.  Perhaps  the  best  graphical  examples  of  this  phenomenon  are  the 


x(t) 


Fig,  2:1:3   Sine  wave  with  amplitude  A  and  period  p 


23 


sine  and  cosine  wave  forms. 

There  are  several  familiar  causes  of  cyclical  variations  in  an 
economic  time  series.   Davis  (3,  26)  mentions  some  of  these,  a  few  of 
which  are  mentioned  below.   Perhaps  most  obvious  cause  for  this  phenome- 
non is  found  in  the  seasons  of  the  year.   A  12-month  period  usually 
results  from  this  seasonal  variation.   A  second  common  pattern  in  economic 
data  is  the  "AO-month  cycle"  which  corresponds  to  that  of  our  national 
economy.   Davis  also  mentions  a  15-20  year  "building  cycle"  shown  by  many 
production  series.   We  cannot  expect  these  cyclical  variations  to  be 
completely  uniform,  since  in  nearly  every  case,  we  are  dealing  with 
human  conduct. 

The  theory  of  Fourier  series  provides  the  basis  for  time  series 
representation  by  trigonometric  models.  Fourier  showed  that  any  reason- 
able periodic  function  of  time  can  be  represented  by  taking  a  sufficient 
number  of  terms  in  the  series 

5(t)  =  A^  +  A^cos[^(t-a^)]  +  A2Cos[^(t-a2)] 

+  .  .  .  +  A^cos[^  (t-aj^)]  +  .  .  . 

where 

A  ,  A  ,  A_,  .  ,  .  A  ,  .  .  .  -  constants, 
o   1   2         n 


Applying  the  identities 


r~2 r 

VAC  +  B'    =•  amplitude 


arctan  -rrr  "  phase  angle 
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we  may  transform  our  equation  for  5(t)  Into  a  series  of  sines  and 
cosines. 


at)  =  A^  +  A£  Sin  ^  +  A-  sin  ^  +  .  .  . 


where 


and 


,  _  ..   2kTit  , 
+  Aj^  Sm  -^  + 


+  B'^  cos  ^  +  B^  Cos  ^  +  .  .  . 


+  Br  Cos  ^  + 


p/k  =  the  periods  of  the  cycles  represented  by  each  sine- 
cosine  pair 


K'   ^V   ^2 K'    ^2' 


constants. 


As  before,   we  could  write   this  model   in  the  form 


e(t)   ■=  x(t)    =     E      a    (t)    f    (t) 
1=1 


=  a-    (t)    f(t) 


where 


a'(t)   = 


^2 


and   f(t) 


Sin  2Tit/P 

Cos  2iit/P 

Sin  4Tit/P 

Cos  4Tit/P 
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Again  note  that  the  coefficients  are  time-dependent  to  allow  for  grad- 
ual changes  in  the  process.   Thus  they  must  be  re-estimated  periodically 
as  additional  data  is  made  available.   By  including  two  terms,  sine  and 
cosine,  for  each  cycle  period,  we  allow  the  origin  of  time  and  the  mean 
to  be  arbitrary;  i.e.,  we  allow  for  a  difference  in  phase  angle.   Obser- 
ving the  process  at  discrete  sampling  intervals,  the  highest  frequency 
that  can  be  detected  is  the  Nyquist  frequency  such  that  there  are  two 
observations  per  cycle.   Brown  (2,  74)  warns  us  that  trigonometric 
models  should  be  used  only  where  there  is  a  known  underlying  cause  for 
a  periodic  time  series  form. 

Having  described  the  basic  mathematical  time  series  forms  in  their 
pure  states,  we  move  now  to  the  consideration  of  more  practical  examples. 
This  brings  us  to  Brown's  third  time  series  model  classification,  com- 
posite models.   These  present  nothing  new,  except  the  combination  of 
the  previous  models,  algebraic  and  transcendental.   The  methods  of 
exponential  smoothing  make  it  possible  to  estimate  coefficients  for  any 
model  composed  of  terms  that  are  sums  or  products  of  algebraic  and 
transcendental  functions  of  time,  in  any  combination  whatsoever. 

An  example  may  best  serve  to  illustrate  the  utility  of  such  a  com- 
posite model.   Consider  the  hypothetical,  though  typical,  curve  shown 
in  Figure  2:1:4. 


26 


x(c) 


Fig.  2:1: A  Hypothetical  time  series  plot 

There  are  several  obvious  elements  in  the  process  which  describes  this 
time  series.   First,  the  model  will  probably  contain  constant  and  linear 
terms  to  show  the  starting  point  and  apparent  linear  trend  of  the  series. 
In  addition  to  these  terms,  a  pair  of  cyclic  terms  with  period  k  will 
probably  be  necessary.   Since  the  amplitude  of  the  cyclic  variation 
appears  to  be  increasing,  we  may  want  to  Include  the  product  of  a  poly- 
nomial term  and  a  pair  of  trigonometric  terms.   Thus  far,  we  have 
formulated  the  following  representation  for  the  time  series,  x(t). 

x(t)   ■=  tCt)   +  e(t) 
=  x(t)  +  e(t) 


=     Z     a   (t)   f.(t) 
i-1 

^_,      2iTt    ,  „        27It 

"  a,    +  a.t  +  a.    Sin  — ; ^■  a,    Cos  — ; — 

12  3  k  4  k 

„,         27lt  „  2TIt 

+  a.t   Sin  -; —  +  a.t  Cos  ^; — 

5  k  6  k 
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2-nt 

a^  +  a^t  +  (^^  +   ajt)  Sin  -^ 
+  (a,  +  a.t)  Cos  —r—  • 

HO  IC 


where  4(t)  =  x(t),  if  our  model  is  the  correct  one 

a.  =  constant  for  i  =  1,  2,  .  ,  .,6,  for  the  short  run. 

In  actual  practice,  a.  =  a  (t) ;  that  is,  we  must  estimate  the  fitting 
function  coefficients  from  our  smoothing  process.   So,  ^ (t)  =  x(t)  only 
if  we  estimate  a,(t)  accurately  and  our  model  is  the  correct  one.   If 
it  is  felt  necessary,  additional  terms  can  be  included  in  the  process 
model,  such  as  a  higher  harmonic  to  the  cyclic  component.   It  may  seem 
at  first  glance  that  we  should  try  to  describe  J (t)  as  accurately  as 
possible;  however,  as  Brown  notes,  for  each  degree  of  freedom  added  to 
the  model  the  work  required  to  estimate  the  coefficients  will  approxi- 
mately double,  and  the  computing  time  will  increase  proportionally  to 

3 
n  . 

The  fourth  and  last  of  Brown's  time  series  model  classifications 

is  regression  models.   This  classification  covers  a  wide  group  of  linear 

forecast  models  of  the  form 

5(t)  =  a,f,  (t)  +  a,f,(t)  +  .  .  .  +  a  f  (t) 
11      z  /  n  n 

where  f  .  (t)  can  be  any  arbitrary  function  at  all.   Many  times  the 
fitting  functions  will  be  one  or  a  combination  of  "leading"  series 
related  to  x(t).   Brown  offers  several  arguments  against  the  use  of 
regression  models.   First,  there  seems  to  be  no  statistical  test  which 
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can  be  applied  to  the  model.   Secondly,  one  usually  doesn't  know  the 
independent  series  long  enough  in  advance  to  be  useful  in  a  forecast.      : 
Finally,  the  computational  techniques  involved  are  considerably  less 
efficient  than  those  of  exponential  smoothing  with  algebraic  and  tran- 
scendental types  of  fitting  functions.   As  we  will  show  in  a  later  section, 
the  efficiency  of  computation  in  exponential  smoothings  depends  on  the 
mathematical  relationships  between  fitting  function  values  at  different 
sampling  intervals.   We  can  not  establish  such  relationships  when  the 
fitting  functions  are  arbitrary  time  series. 

2.2  Trend  Analysis 

By  trend,  or  secular  trend  as  it  is  more  commonly  called,  we  mean 
"...that  characteristic  of  the  series  which  tends  to  extend  consistently 
throughout  the  entire  period"  (Davis,  3,  15).   The  period  referred  to  in 
this  definition  is  the  total  time  covered  by  the  available  data.   Trend 
analysis  is  an  attempt  to  determine  this  secular  trend  in  the  time  series. 
It  is  performed  on  the  data  for  several  reasons: 

1.  To  determine  the  polynomial  component  of  the  time  series  model. 

2.  To  aid  in  detrending  the  time  series  for  the  auto  correlation 
and  spectral  analyses  to  be  discussed  in  the  next  two  sections 
of  this  chapter.   By  detrending  we  mean  to  reduce  the  series 
in  such  a  manner  that  the  residual  fitting  errors  will  have  an 
expected  value  of  zero. 

3.  To  aid  in  identification  of  any  amplitude  growth  that  may  be 
occurring  within  the  cyclic  components  of  the  time  series. 
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Davis  goes  on  to  describe  four  of  the  more  commonly  encountered 
trend  forms  in  economic  time  series  analysis.   The  first  of  these  is  the 
straight  line  trend,  described  by  a  first  degree  polynomial  on  t,  time. 
Another  common  trend  type  is  a  straight  line  fitted  to  data  plotted  on 
a  logarithmic  scale.   This  type  corresponds  to  the  simple  exponential 
model  discussed  earlier.   He  calls  the  third  trend  line  type  a  "logistic 
curve"  or  "curve  of  growth".   This  curve  is  more  applicable  to  a  descrip- 
tion of  the  growth  of  new  industries  and  population  series.   It  may  be 
described  as  a  transition  line  between  a  lower  initial  level  of  x(t),  to 
a  stable  "saturation"  point  on  a  higher  level.   The  logistic  curve  is  repre- 
sented mathematically  by 

Ut)  =  - 


where 


1  +  e-h' 


K  ,  K,  =  constants 

C(t)  =  only  the  trend  component  of  the  time  series  in  this 


The  form  of  the  logistic  curve  is  shown  in  Figure  2:2:1. 


Ut) 


Fig.  2:2:1  Example  of  logistic  curve 
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The  last  of  the  trend  lines  described  by  Davis  is  the  familiar 
moving  average . 

We  may  consider  C  Ct)  as  being  composed  of  two  components,  trend 
and  periodic. 

C(t)  =  trend  +  periodic 
=  T  +  P 
where 

T  =  the  trend  component 
P  =  the  periodic  component 

and  both  T  and  P  are  deterministic  functions  of  time.  We  wish  to  describe 
T  so  that  we  may  study  P  more  closely  by  performing  the  following  operation: 

x(t)  =  e(t)  +  e(t) 

=  T  +  P  +  e(t) 
P  +  c(t)  =  x(t)  -  T 

As  noted  earlier,  the  periodic  component  is  expressed  as  a  series  of  sines 
and  cosines.  These  trigonometric  functions,  however,  may  be  expanded  into 
a  series  of  polynomial  terms. 


Sin  y  =  y  -  ^  +  ^ 


Zl  +  2_  _  .  .  . 


where 


cosy  =  l-^  +  .^j 


2k¥t     ,  , 

y  = as  before. 

P 

Here  we  can  see  the  danger  of  mixing  the  trend  and  periodic  components, 
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especially  when  high  degree  polynomials  are  used  to  describe  trend.   Unless 
we  are  cautious,  we  may  remove  part  of  the  periodic  variation  when  we  are 
detrending  the  time  series.   This  would  probably  lead  to  grave  forecast 
errors  later  on.  A  plot  of  the  trend  line  with  the  actual  observations 
may  aid  the  investigator  in  avoiding  this  pitfall. 

In  order  to  illustrate  the  detrending  process,  let  us  consider  an 
ideal  model  of  the  form 


271 1 


x(t)    =  5(t)    =  a,    +  a.t  +    (a,  +  a^t)    sin 

L  Z  Jo  p 


.      /  ,  N  2TVt 

+    (a,    +  a,t)    cos  — — 
4  6  p 


where  £(t)  is  always  zero.  This  time  series  is  shown  graphically  in 
Figure  2:2:2. 


<(t) 


Fig.  2:2:2  Ideal  time  series  with  no  random  variation 


In  this  example  we  may  identify  the  process  components  as 
T  =  a^  +  a^t 


P  =  (a„  +  a^t)  sin  +  (a.  +  a.t)  cos  

J    J        p      4    o        p 
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The  detrended  data  will  have  the  form  of  P  alone,  as  shown  in  Figure  2:2:3. 


Fig.  2:2:3  Detrended  form  of  ideal  time  series  above 

The  series  now  has  an  expected  value  of  zero.   Thus,  it  is  easier  to 

Identify  the  cycles  present  in  the  data  and  any  growth  in  amplitude  that 
may  be  occurring.   The  detrended  data  may  then  be  subjected  to  autocor- 
relation and  spectral  analysis  to  determine  the  periodic  components 
present.   Polynomial  regression  (Appendix  C)  is  perhaps  the  best  method 
to  use  in  detrending  the  data,  and  it  is  the  method  used  in  this  research. 


2.3  Autocorrelation  Analysis 

Though  "detrending"  aids  in  the  identification  of  periodic  compon- 
ents, it  is  usually  still  inadequate.   Random  noise  and  the  superposition 
of  one  harmonic  upon  another  often  make  it  impossible  to  recogniEe  all 
of  the  periodic  components  present  in  the  process  involved.   Autocorrel- 
ation analysis  is  one  method  of  manipulating  the  data  to  indicate 
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(1)  Wliether  or  not  a  process  exists  in  the  detrended  data 

(2)  Whether  or  not  the  data  contains  periodic  components 

(3)  At  least  the  largest  periodic  component  detectable  in  the  data 

(4)  Any  Icig  relationship  between  data  values  a  certain  time  period 
apart  in  the  time  series 

By  applying  the  "detrending"  techniques  of  the  previous  section, 
our  data  is  now  in  the  form 

x^(t)  =  P  +  e(t)  =  x(t)  -  T 

where  x.(t)  =  the  detrended  form  of  the  data  and  P  and  T  are  the  periodic 
d 

and  trend  components  of  the  time  series,  respectively.   Since  x  (t)  now 
has  an  expected  value  of  zero,  by  the  principal  theorem  of  Fourier  we 
may  write 

x^(t  +  k)  ;  x^(t) 

where  k  =  the  period  of  the  fundamental  periodic  component  of  the  series. 

This  equation  assumes,  of  course,  that  a  periodic  component  exists, 

without  any  amplitude  growth.   The  fundamental  periodic  component  here  is 

the  largest  detectable  period  in  the  available  time  series  data. 

Autocorrelation  analysis  utilizes  a  function  known  as  the  average 

lagged  product,  which  may  be  written  as 

T 

t=k+l  ^   '^  '^ 
T  -  1  -  k 

where  k  is  the  lag,  or  period  of  time  between  the  time  series  values 
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whose  products  we  are  summing.   This  average  lagged  product,  performed 
on  a  sequence  of  numbers  whose  expected  value  is  zero,  is  the  autocovar- 
iance  for  lag  k,  R   (k) . 


I   X4l(t)x,(t-k) 

XX 


T  -  1  -  k 

The  variance  of  a  sequence  of  numbers  whose  expected  value  is  zero  is 

T 
ox/  =  ^1^  x^(t)VN 

where  N  =  the  total  number  of  observations. 

This  is  just  Rxx(O),  the  average  lagged  product  for  a  lag  of  k  =  0. 

2 
Rxx(O)  =  0  . 
xd 

Note  here  that  the  function  Exx(k)  will  always  be  symmetrical,  so  that 

Rxx(-k)  =  Rxx(k) 

The  normalized  form  of  the  autocovariance  function  Rxx(k)  is  simply 
„,,  ,  _  Rxx(k) 

Since  Rxx(O)  >  Rxx(k)  for  any  k, 

-  1  <   (k)  <  +  1 
and 

P(0)  =  +  1. 
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For  essentially  random  noise  p (k)  will  be  nearly  zero  for  all  k^O,   and 
p(0)  =  1.   So  if  no  periodic  component  exists  in  the  data,  a  plot  of 
the  autocorrelation  function  against  the  lag  k  would  appear  as  shown  in 
Figure  2:3:1. 


Fig.  2:3:1  Autocorrelation  function  for  time  series  with  no 
period  component 

p(k)is  shown  in  Figure  2:3:1  as  p (k)  to  indicate  that  this  is  a  computed 
value  subject  to  random  variation.   Since  Rxx(k)  is  symmetrical,  P (k) 
will  also  be  symmetrical,  and  we  need  only  plot  P (k)  for  k>0  to  study 

the  autocorrelation  function.   The  standard  deviation  of  the  computed 

-I' 
autocorrelation  coefficients,  P (k) ,  is  a  =  (T-k-l)'^  for  a  sample  of  T 

observations.   Brown  warns  us  that  the  autocorrelation  coefficients 

should  be  treated  as  zero  unless  they  are  clearly  greater  in  magnitude 

than  the  standard  deviation,  a  .   He  also  says  that  the  maximum  lag  should 

not  exceed  kmax  =  T/20.   This  restriction  may  be  too  tight  to  follow  in 

practice,  since  to  study  a  yearly  cycle  this  means  that  at  least  twenty 

years  of  monthly  data  would  be  necessary. 

It  may  be  helpful  at  this  point  to  discuss  the  possible  interpreta- 
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tions  that  may  be  made  from  the  autocorrelation  function  plot.   Consider 
the  plot  shown  in  Figure  2:3:2.   The  pattern  shows  a 


1.0 


P(k)    o 


-1.0 


Fig.  2:3:2  Example  of  an  autocorrelation  function 
trough  at  six  months,  a  peak  at  twelve  months,  and  inflection  points 
at  three  months  and  nine  months.   This  is  indicative  of  the  existence 
of  a  twelve  month  periodic  component  in  the  data.   The  squared  nature 
of  the  trough  and  crest  are  indications  that  higher  harmonics  may  be 
present  also.  Consider  now  the  plot  in  Figure  2:3:3.  Notice  again  the 
apparent  twelve  month  periodic. 


1.0 


p(k)    0 


-1.0 


12       15 


Fig.  2:3:3  Example  of  an  autocorrelation  function 
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However,  notice  also  the  smaller  peak  at  six  months  and  the  cosine 
wave-like  form  of  the  curve  between  k  =  1  and  k  =  6.   This  Indicates 
that  a  six  month  harmonic  may  also  be  included  in  the  periodic  component 
of  the  process. 

We  may  sunmiarize  the  autocorrelation  analysis  procedure  in  the 
following  steps: 

(1)  Detrend  the  time  series. 

(2)  Compute  autocorrelation  function,  P (k) ,  for  the  desired  lag 
values,  k. 

(3)  Plot  k  versus  P (k) . 

(4)  Analyze  and  interpret  the  plot  in  (3) . 

It  may  be  noted  here  that  the  more  completely  we  detrend  the  data, 
the  sharper  will  be  the  peaks  on  the  autocorrelation  function  plot,  and 
hence,  the  more  easily  we  may  recognize  the  existing  periodics.   We  must, 
of  course,  be  careful  not  to  remove  any  periodic  variation  by  mistaking 
it  for  trend. 

Autocorrelation  analysis  is  especially  advantageous  when  the  time 
series  may  be  described  using  fitting  functions  which  are  previous  values 
of  the  same  time  series  k  observations  in  the  past.   This  is  referred  to 
as  an  autoregressive  model.   For  time  dependent  periodic  components 
another  analytical  technique  has  attributes  not  offered  by  autocorrelation 
analysis.   This  latter  technique,  harmonic  analysis,  is  discussed  in  the 
following  section. 

2.4  Harmonic  Analysis 

Recall  from  Section  2.1  that,  according  to  the  fundamental  theorem 
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of  Fourier  series,  any  reasonable  periodic  function  of  time  with  zero 
mean  can  be  represented  arbitrarily  closely  (Brown,  2,  396)  by  taking 
terms  from  the  infinite  Fourier  series 

,(t)=|A^+  I     A^Cos(Sf)+  "  B^SlnC^). 
n=l  n=l 

-  a  <  t  <  a 


where 


A  ,  B  =  constant  coefficients 
n   n 


A  =  constant  (zero  for  the  detrended  data) 
o 

a  =  one  half  the  range  of  the  main,  or  longest,  period, 


The  coefficients  A  and  B  may  be  computed  from  the  following  integrations: 
n      n 

A  =  i  /  x(t)Cos(22^)dt 
n   a  a 


B  =  i  ;  x(t)Sin(J^)dt 
n   a  a 


The  period,  p,  and  amplitude  R,  of  the  nth  harmonic  term  may  be  written 


2a 

?„  =  -;: 


=  v/?T 


The  periods  represented  in  this  Fourier  series  make  up  the  Fourier 
sequence ,  Chat  is , 
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2a,  2a,  2a 2a, 

12   3         n 


The  advantage  of  the  Fourier  sequence  is  revealed  in  the  following 
identity: 

R^  +  R^  +  R3  +  •  •  •  +  ^n  +  ■  •  •  =  2" 

where  o   =  the  variance  of  the  data  series,  5 (t) . 

Thus  we  see  that  the  ratio  R^/2a   gives  us  the  fractional  variation, 
or  energy,  explained  by  the  inclusion  of  the  nth  harmonic  term  in  our 
time  series  model.   Our  goal  in  describing  the  periodic  components  of 
the  data  is  not,  however,  to  fit  a  Fourier  series  to  the  data.   Rather 
we  wish  to  determine  the  periods  which  the  time  series  possesses  and, 
if  possible,  determine  the  contribution  of  each  of  these  periods  in 
describing  the  time  series.   It  is  conceivable  that  one  or  more  of  these 
periods  would  not  be  included  in  the  Fourier  sequence.   Hence  we  may 
define  harmonic  analysis  as  the  technique  of  determining  the  principal 
harmonic  elements  of  a  given  set  of  data. 

For  our  purposes,  a  better  sequence  of  periods  to  study  is  the 
arithmetic  sequence 

1,  2,  3,  4 a. 

This  sequence  may  also  be  used  to  describe  the  time  series  data,  though 
the  energies  for  these  periods  are  no  longer  additive.   We  may  now 
describe  the  time  series  in  the  form 


40 


£(t)=|A  +  t     ApCosC^)  +  f  BpSln(^) 
p=l        ^     P=l 

where  A  ,  B  =  the  coefficients  of  the  harmonic  terras  for  a  period  of  p. 
P   P 

It  is  helpful,  for  reasons  to  be  discussed  later  in  this  section,  to 

subject  the  detrended  data  to  the  harmonic  analysis  rather  than  the 

raw  time  series  data.   Since  the  detrended  data  may  be  written 

x,(t)  =  P  +  e(t) 
d 

and  the  periodic  component,  P,  has  a  zero  expected  value,  we  may  now 

write 

P  =  E  A  Cos(^  +  I      B  Sin(^). 
P=l  P     P     p=l  P     P     . 

We  may  determine  the  amplitudes  of  the  periodic  terms  from  the  following 
formulas : 

A  =  |t  I  x,(t)Cos(^) 

P   ^  t=l  P 

■>       ^'  2  t 

B  =^  E  x,(t)Sln(-— ) 
p   N  j.^^  d        p 


\=  / 


A^TT^ 
p     p 


where  N'  is  chosen  equal  to  the  largest  integral  multiple  of  p  In  the 
total  number  of  observed  data  values,  N. 

Having  determined  the  values,  R,  of  the  amplitudes  of  the  harmonics 
in  the  arithmetic  sequence,  we  wish  to  compare  their  magnitudes.   This 
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will  give  us  a  relative  measure  of  the  contribution  of  each  harmonic 
to  the  description  of  x^(t).   The  most  widely  used  method  of  compaisson 

is  the  "periodogram,"  introduced  by  Sir  Arthur  Schuster  (1851-193A) 

2 
(Davis,  3,  31).   A  "periodogram"  Is  the  graph  of  R^  or  R^  against  p. 

Graphs  of  R  versus  p  have  been  employed  throughout  this  research. 

This  graph  will  contain  peaks  above  the  periods  corresponding  to 
harmonics  present  in  the  data,  x^(t).   Generally,  the  higher  the  peak, 
the  more  significant  is  that  particular  harmonic  in  describing  the 
periodic  component  of  the  process.   Hence,  inclusion  of  that  harmonic 
in  the  time  series  model  will  account  for  a  greater  decrease  in  error 
variation  than  would  a  harmonic  with  a  lower  peak. 

Consider  the  hypothetical  "periodogram"  shown  in  Figure  2:4:1. 
From  this  graph 


300 


R    200 
P 


100 


J p  (months) 


12     4     6     8     10    12 
Fig.  2:4:1  Hypothetical  harmonic  analysis  periodogram 
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we  may  infer  that  a  periodic  component  with  P  =  12  months  should  probably 
be  included  in  our  time  series  model.   Also,  a  six-month  harmonic  may  be 

desired  in  the  model. 

We  may  obtain  some  idea  of  the  effect  of  the  inclusion  of  certain 

periodic  terms  in  the  time  series  model  by  computing  a  per  cent  reduction 
in  error  variance  for  x^(t).   For  example,  let  a^,  be  the  variance  of 
the  detrended  data,  and  let  0^.,  be  the  variance  of  detrended  data  minus 
the  periodic  terms  we  wish  to  include.  For  the  example  in  Figure  2:4:1 
this  would  be:  . 

x,(t)  =  p  +  E(t) 

d  -  ■■  :   ,  .1  < 

x/(t)  =  x^(t)  -  P' 


where 


p'  =  Ag  cos  (^)  +  B^  sin(-g-) 

+  A^2  Cos  (^)  +  Bi2  ""^TI^ 


Now  we  may  write 

2 

°x   d 
E  ,=  100[1  -  H5 — ] 

xd 


where 


E   =  the  per  cent  reduction  in  0  ,  due  to  P'. 
P' 


Ideally  we  would  like  our  "periodograms"  to  exhibit  sharp  peaks 
for  the  harmonics  that  should  be  included  in  the  time  series  model. 
Since  our  data  still  contains  some  random  variations,  sharp  peaks  are  not 
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always  possible.   By  working  only  with  the  detrended  data,  we  lessen 
the  chances  of  misinterpretation  of  the  graph.   Inclusion  of  the  trend 
component  in  the  harmonic  analysis  may  give  erroneous  peaks  on  the 
graph,  or  hide  harmonics  that  should  have  been  Included  in  the  model. 
Even  after  trend  is  removed  from  the  data,  careful  interpretation  of 
the  "periodogram"  is  essential.   A  very  broad  peak,  such  as  that  shown 
in  Figure  2:4:2,  suggests  the  existence  of  more  than  one  harmonic. 


300 


200 


100 


16  12 

Fig.  2:4:2  Example  of  periodogram  with  broad  peak 

Existing  periods  of  very  low  frequency  may  often  be  handled  adequately 
by  the  adaptive  nature  of  the  smoothing  process. 

We  may  see  from  this  discussion  that  harmonic  analysis  not  only 
identifies  the  harmonics  existing  in  the  time  series,  but  also  indicates 
their  relative  contributions  in  describint  the  process.   Nevertheless, 
it  is  still  subject  to  error  and  misinterpretation.   Probably  a  wise 
procedure  is  to  perform  both  the  autocorrelation  and  harmonic  analyses, 


AA 


so  that  one  may  serve  as  a  check  on  the  other. 


■  -,i^.— -.n--:.-^-J 
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3.   FORECASTING  THE  TIME  SERIES 

3.1  Justification  for  Exponential  Smoothing 

We  have  derived  a  simplified  formulation  of  the  exponential  smoothing 
process  in  section  1.3  of  Chapter  1.   In  deriving  the  expressions  shown 
there,  we  were  dealing  with  the  simplest  of  time  series  models,  the 
constant  process.   It  was  shown  that  the  function  S^(x)  was  a  good  estimate 
of  the  value. of  the  process  at  time  (t+1) .   For  this  model  we  had  only  one 
fitting  function,  t^it)   =  1,  where 
Kit)    =  a(t)  f^(t) 

=  a(t)  (1) 
x(t)  =  a(t) 

Thus  we  see  that  once  we  have  formulated  a  good  model  of  the  true  process, 
we  use  exponential  smoothing  to  estimate  the  coefficients,  a^(t) ,  of  the 
fitting  functions  in  the  model. 

To  bring  us  up  to  date,  let  us  briefly  review  our  progress  at  this 

point: 

(1)   The  time  series  may  be  described  as 

x(t)  =  T  +  P  +  e(t) 


where 


T  =  the  trend  component 
P  =  the  periodic  component 

e(t)  =  the  random  variation.  Both  T  and  P  are  determin- 
istic functions  of  time. 
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(2)  By  means  of  polynomial  regression  we  may  determine  the  fitting 
functions  in  T. 

(3)  We  may  determine  the  fitting  functions  in  P  through  autocorrela- 
tion and  harmonic  analyses. 

(4)  Now,  through  exponential  smoothing,  we  may  estimate  the  coeffi- 
cients of  all  the  fitting  functions  in  the  time  series  model, 
and  then  make  forecasts  of  future  values  of  the  time  series. 

From  our  previous  development  of  exponential  smoothing  we  arrived 
at  the  expressions 

S^M    =  aXj,  +  BS^_^(x) 

2 
=  «Xj.  +  Ci3x^_^  +  aB  Xj._2  +  .  .  . 

+  ag"x    +  .  .  .  +  e*^  X 
t-n  o 


S^(x)  =  a\  S'^  x^_^  +  6\ 
n=0 


where 


S  (x)  =  the  exponentially  smoothed  value  of  x(t)  up  to 

time  t. 
a  =  smoothing  constant 
S  =  1-a  =  discount  factor. 

Note  in  this  expression  once  again  the  advantages  offered  by  exponential 
smoothing: 

(1)   Data  is  weighted  to  reflect  decreasing  importance  with  increasing 

age. 
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(2)  All  past  data  is  included  in  the  statistic  S  (x) ,  where  t  is 
taken  to  be  the  present  time. 

(3)  Storage  of  past  data  requires  only  the  retention  of  the  current 
value  of  S  (x) . 

As  a  criterion  for  justification  of  the  exponential  smoothing  process, 
we  shall  say  that,  if  the  expected  value  of  S(x)  equals  the  expected 
value  of  X,  the  procedure  is  justified.   Thus,  S  (x)  is  a  valid  method 
for  estimating  a. (t) ,  the  constant  process  fitting  function  coefficient. 
So,  we  may  say 

E[S(x)]  =  E[«  I   b'^x^  ,  ]  +  E[e"x  ]   . 
k=0    """         ° 


The  last  term  shown,  EfS^x  ],  is  zero,  since  S  is  a  number  less  than  one. 
o 


Now, 


E[S(x)]  =  E[a  I      b\        ] 
k=-0    '^"'' 


a  Z   b''  E[x   ] 
k=0       ^'^ 


E[x]a   E   6^ 
k=0 


Z   B  is  the  summation  of  a  converging  geometric  series.   It  has  the 
k=0    ^ 
value  rr^     for  an  infinite  number  of  terms.   So, 

E[S(x)]  =  E[x]  a  Jig 
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=  E[x]  - 

a 

E[S(x)]  =  E[x] 

and  the  procedure  is  justified.   In  addition  to  the  advantages  cited 
above,  the  response  rate  of  exponential  smoothing  is  quite  flexible  and 
the  computations  simple,  as  will  be  shown  in  later  discussions.   Note  also 
another  characteristic  of  the  smoothing  process.   Let  us  rearrange  the 
previous  expression  for  S  (x)  as  shown  below. 

Sj.(x)  =  ax^  +  3S|._^(x)  ... 

=  aXj,  +  (l-a)Sj._^(x) 

=  Sj._^(x)  +  a(X|.  -  S^_^(x)) 

The  second  expression  on  the  right  is  merely  the  current  observation,  x  , 
minus  the  estimate  of  x  at  time  t-1,  S   -,  (x) ,  all  multiplied  by  the 
smoothing  constant.   So,  this  is  a  measure  of  the  error  of  the  previous 
forecast.   Hence,  the  smoothing  procedure  updates  the  coefficient  esti- 
mates to  adjust  for  the  previous  estimation  error.   The  accuracy  of  the 
procedure  may  be  inferred  from  the  fact  that  the  function  S(x)  minimizes 
the  weighted  sum  of  squared  residuals  (Appendix  B) .       ■   ■ 

3.2.   Discounted  Multiple  Regression 

Until  now,  we  have  been  working  with  a  time  series  that  may  be  repre- 
sented by  a  constant  process.  We  now  wish  to  consider  a  more  practical 
example,  where  the  process  is  composed  of  an  arbitrary  set  of  functions. 
For  these  fitting  functions  we  determine  values  of  the  coefficients  that 
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miniciize  the  weighted  sum  of  square  residuals,  as  before. 

Consider  again  the  general  process  mentioned  in  section  1.3. 


S(t)  =  a^(t)f^(t)  +  a^Mi^it)  +  . 


+  a  (t)f  (t) 
n    n 


Z      a  (t)f(t) 

i=l 


where 


£(t)  =  the  process,  or  the  time  series  minus  random  noise 
f . (t)  =  the  i   fitting  function  evaluated  at  time  t 
a.(T)  =  the  i   fitting  function  coefficient  estimated  at 
time  T. 
The  fitting  functions  may  be  any  deterministic  functions  of  time,  as 
long  as  their  values  are  known  exactly  for  the  present  and  for  the  time 
in  the  future  for  which  a  forecast  is  desired.   This  restriction  allows 
the  inclusion  of: 

(1)  Simple  mathematical  functions  of  time. 

(2)  Empirical  functions  that  satisfy  the  restriction. 

(3)  Previous  observations  of  the  dependent  time  series. 

Having  determined  the  fitting  functions  to  be  used,  we  must  estimate  their 
coefficients  so  that  the  desired  x-value  may  be  obtained  from 

x(t)  =  .^^  a.(t)f .(t) 
=  a'(t)f(t) 
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where  a(t)  and  f(t)  are  n-component  column  vectors  evaluated  at  time  t. 
Recall  that  the  process  C (t)  may  change  slightly  over  a  long  period  of 
time  without  injecting  significant  error  into  the  smoothing  process. 
Over  the  short  run,  however,  we  may  treat  the  coefficients  as  constants 
and  obtain  a  least  squares  estimate  of  their  values  from  the  exponentially 
discounted  data. 

The  forecast  of  x  for  ^    periods  into  the  future,  considering  the 
present  time  as  t,  would  then  be 


x(t+T)  =  a^(t)f^(t-K)  +  a^Mt^(.t-t^) 
+   .    .    .  +  a  (t)f„(t+T) 


i!l  a.(t)fi(t+^) 


a'(t)f(t+T) 


where 

a.(t)  =  a  least  squares  estimate  of  a.(t)  based  on  the 
weighted  data,  x,  from  t  =  0  up  to  the  present 
time,  t. 
f  (t+x)  =  the  i'   fitting  function  evaluated  at  the  future 
time,  t+''^  for  which  the  forecast  is  desired. 
We  are  assuming  that  the  estimates,  a  (t) ,  based  on  all  past  data,  are  the 
best  that  we  can  obtain.   The  least  squares  estimate  of  a^(t)  is  similar 
to  a  multiple  regression  solution;  however,  in  this  case,  we  are  minimizing 
the  discounted  sum  of  squared  residuals 
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I      3-'[x(T-j)  -  E  a.(T)f  (T-j)] 


where 

T  =  present  time. 
Hence,  we  call  this  procedure  "discounted  multiple  regression. 

The  ensuing  discussion  will  show  the  technique  for  finding  these 
least  square  estimates  of  the  coefficients.   First,  the  following  defin- 
itions should  aid  the  understanding  of  the  procedures  that  follow. 

X  =  [x^,   x^,  Xj x^] 

=  all  past  observed  time  series  values  up  to  the  present 

time  T. 
=  al  X  T  row  vector. 


f(t) 


fj_(t) 
f2(t) 


f^(t) 


an  n  X  1  column  vector  of  the  n  fitting 
functions  evaluated  at  time  t. 


a(t) 


a^(t) 
a^Ct) 


.(t) 


an  n  X  1  column  vector  of  fitting 
coefficients  for  time  t. 
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5^(t) 


fj^(l)    f^(2) 


f^m 


f^d)    £,(2)  .  .  .  f2CT) 


f^(l)    f^(2)  .  .  .  f^(T) 


=  an  X  T  matrix  with  elements  f.(t)  in  the  i   row  and  the 

th 
t   column. 

a(t)  is  our  least  squares  estimate  of  a(t)  and  is  of  the  same  vector  form 
as  that  shown  above  for  a(t).   Since  we  have  defined  the  time  series 
estimate  as 

x(t)  =  a'(t)  f(t) 

we  may  denote  the  corresponding  residual  as 

e(t)  =  x(t)  -  x(t) 

=  x(t)  -  a'(t)f  (t) 

Now  the  discounted  sum  of  squared  residuals  would  be 

T 

Z     w^t)e^t) 
t=l 


where 


w(t)  is  the  t  diagonal  element  of  a  T  x  T  matrix  W(t) . 


W(t)  = 


w(l) 


w(2) 


w(T) 
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where 


w^  (T-j)  =  B^  j  =  0,  1,  2,  .  .  .,  T-1. 


Now,  it  we  define  a  row  vector  of  residuals,  e(t), 
e(t)  =  x(t)  -  a'(t)S^(t), 

then  the  discounted  sum  of  squared  residuals  Is 

S  =  [e(t)W(t)]  [e(t)W(t)]' 
a 

=  [x(t)W(t)  -a'(t)^(t)W(t)]  [x(t)W(t)  -a'(t)5'(t)W(t)]' 

for  a  particular  set  of  coefficients  a(t).  To  minimize  this  sum,  we  need 
to  select  a(t)  such  that 


^^^  =  0. 


3a(t) 

This  expression  gives  us  n  simultaneous  linear  equations  which  can  be 
solved  for  the  n  coefficients  a.(t).   As  shown  in  Appendix  B,  the  solution 
to  this  expression  is 

a'(t)  =  x(t)W(t)W'(t)?''(t)F"^(t) 
where 

F(t)  =^(t)W(t)W'(t)?''(t) 

=  an  n  X  n  syiraaetrical  matrix 

As  long  as  T>>n,  F(t)  will  have  an  inverse.   If  we  define  a  data  vector 
g(t)  as 

t'(t)  =  x(t)W(t)W'(t)5^(t) 
then 

a'(t)  =  g'(t)F"''-(t). 
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Note  that  F(t)  does  not  depend  on  the  data.   It  does,  however,  depend  on 
time,  and  must  be  updated  at  each  sampling  interval.   Looking  at  one 
element  in  the  matrix  F(t) , 

T-1 

Z 

where 


i-i 
Fik(T)  =  I      B^f^(T-j)fj^(T-j) 


T  =  present  time,  we  see  that  It  may  be  updated  by 

F.k(T+l)  =  f  .(T+l)f,  (T+1)  + 

T-1 
I     B^f  .(T-j)f,  (T-j) 

So,  generally,  F(t)f'(t)  +  BF(t-l) 

Since  F(t)  5^  F(t-l),  then  F~-'-(t)  ?!  F"-'-(t-l)  and  a  new  f"  must  be  computed 

at  each  sampling  interval.   Reference  to  an  expression  for  the  i   data 

vector  component  at  time  T, 

T-1  , 

I 
j=0 


T-1   . 
g.(T)  =   I   S^x(T-j)f  (T-j), 
1        j_n 


indicates  that  g(t)  may  be  updated  at  each  sampling  interval  by 

g(t)  =  x(t)f  (t)  +  6g(t-l) 
where 

t  =  the  new  present  time  value 

t-1  =  the  previous  time  period. 
So,  having  updated  F(t)  and  g(t)  and  after  a  sufficient  amount  of  data 
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has  been  gathered  to  allow  the  existence  of  F~  (t) ,  we  may  estimate  a 


by 


a'(t)  =  gCt)F"^(t) 


and  make  a  forecast  of  future  observations  by 

xCt+T)  =  a-(t)f(t+T). 

3.3  Construction  of  the  Transition  Matrix 

We  can  simplify  the  computations  of  the  previous  section  considerably 
If  we  place  a  restriction  on  the  class  of  fitting  functions;  that  is, 
the  vector  of  fitting  function  values  at  time  t+1  must  be  a  linear  combin- 
ation of  the  values  of  these  functions  at  time  t.   Specifically,  we  can 
define  an  n  x  n  matrix  L  with  components  L   such  that 

f^(t  +  1)  =  L^^f^(t)  +  L^2^2^t)  +  .  .  .  +  L^^f  Jt) 
f^Ct  +  1)  =  L^.f^Ct)  +  L^^f^Ct)  +  .  .  .  +  L^^f^Ct) 


^n^^  +  "  =  ^1^1^'^  ■"   "n2^2(^)  +  •  •  •  +  \,f,M 


f(t  +  1)  =  L  f  (t) 

The  transition  matrix,  L,  must  have  an  inverse,  as  will  be  shown  later. 
The  only  fitting  functions  which  satisfy  this  restriction  are 

(1)  Polynomials 

(2)  Exponentials 

(3)  Sinusoids 


56 


once  we  know  £ (o)  and  L.  we  may  find  f(t)  for  any  t  since, 

f(t)  =  Lf(t-l) 

=  L[Lf(t-2)] 

=  L^f (o) 

The  matrix  L  may  be  constructed,  by  inspection  of  the  fitting  functions 
desired,  according  to  rules  regarding  each  allowable  class  of  fitting 

functions . 

For  a  polynomial  of  degree  n,  L  is  an  (n  +  1)  x  (n  +  1)  matrix  with 
ones  on  the  diagonal,  ones  in  the  first  element  to  the  left  of  the 
diagonal  and  zeros  everywhere  else. 

1  0 

1   1 

1   1 


;  ij 


where  L  has  n  +  1  rows  and  n  +  1  columns.   The  mathematics  involved  in 
updating  L  requires  that  polynomial  fitting  functions  now  be  of  the  form 


f  (t)  =  -^  t(t-l)  (t-2) 
n      nl 


. (t-n  +  1) . 


As  mentioned  earlier,  for  a  trigonometric  fitting  function  the  sine 
and  cosine  term  of  each  harmonic  must  both  be  Included  to  account  for  any 
phase  angle  shift.   There  are  two  basic  forms  which  the  corresponding 
transition  matrix  may  take.   For  example,  suppose  we  want  L  for  the 
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trigonometric  pair 


f  (t)  =  sinut 
f^Ct)  =  cos'^t 
L  would  be  the  2x2  matrix 


_^.cosu  sinuj. 

-sintu  costi) 

If  we  want  to  include  this  same  harmonic,  but  with  linearly  growing  ampli- 
tude as  well,  the  fitting  functions  would  be: 

f  (t)  =  sinuit 

f  (t)  =  cosut 

f  (t)  =  t  sinut  ^'  .  '  ■  • 

f ,(t)  =  t  cosut 


and  L  would  be  the  4x4  matrix 


COSOJt 

sinwt 

0 

sinut 

COSLJt 

cosuit 

sinut 

costot 

sinut 

slnut 

cosut 

-sintot 

coswt 

The  elements  of  L  always  take  the  same  position  relative  to  the  diagonal 
of  L,  no  matter  what  the  size  of  L.   So,  for  models  which  include  many 
fitting  functions  or  combinations  of  fitting  functions,  the  transition 
matrix  may  be  constructed  by  stacking  the  simple  matrices  shown  along 
the  diagonal  of  L.   All  other  elements  will  be  zeros.   For  the  construc- 
tion of  exponential  transition  matrices,  refer  to  Brown  (2,  166). 
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3.4  Adaptive  Smoothing 

Observing  the  restriction  placed  on  f (t)  in  the  previous  section 
and  moving  the  origin  of  time  to  the  most  recent  observation,  we  may  now 
write 

xCt+t)  =  £-(T)2(t)  ,  .,  ,-.  ,..   .  '  '  ;. 

This  shift  in  time  is  reflected  only  in  the  fitting  functions,  so  that 
their  values  depend  only  upon  the  increment  of  time  into  the  future  for 
which  a  forecast  is  desired,   x  still  shows  the  actual  time  to  which  the 
forecast  applies,  and  a  shows  the  total  data  upon  which  it  is  based,  all 
data  up  to  time  t.   So,  if  we  shall  always  be  forecasting  one  time  period 
into  the  future,  t  =  1,  we  need  to  compute  f(T)  =  f(l)  only  once. 
The  sum  of  squares  to  be  minimized  now  is 

T 

Z  S^[x(T-j)  -f'(-j)a(T)]^. 

Differentiating  the  weighted  sum  of  squared  errors  with  respect  to  a(t) 
and  setting  the  result  equal  to  zero  as  before,  we  obtain 

I(t)  =  F"^(t)g(t), 


where 


and 


J=0 


t 
g,  (t)  =  Z   B^f  .  (-j)x(t-j) 
j=0 


Again  this  assumes  that  sufficient  observations  have  been  made  for  F(t) 
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to  have  an  inverse,  that  is,  at  least  n  observations  where  F(t)  is  an 

n  X  n  matrix. 

At  each  sampling  interval,  we  must  again  update  g(t)  and  F(t) .   Since 

both  are  functions  of  f(t),  and  f(t)  has  a  different  time  origin  than  in 

discounted  multiple  regression,  the  updating  procedure  will  not  be  the 

same  as  before.   The  data  vector  may  be  expressed 

t 
g(t)  =  x(t)f(o)  +  l      B^f(-j)x(t-j) 

Since 

f(-j+l)  =  Lf(-j),  or  f(-j)  =  L"^f(-j+l), 

t   .   ^ 
g(t)  =  x(t)f(o)  +  Z  B^l"  f(-j+l)x(t-j). 
j  =  l 

If  we  change  the  index  of  summation  to  k  =  j  -  1,  the  data  vector  recur- 
sive relationship  becomes 

g(t)  =  x(t)f(o)  +  BL"-'-g(t-l) 

where  the  coefficients  of  x(t)  and  g(t-l)  no  longer  depend  on  time  and 
thus  become  program  constants.   The  matrix  of  fitting  functions  may  also 

be  computed  recursively  by 

t 

F(t)  =  ::  e^f(-j)f'(-j) 
j=o 

=  F(t-l)  +  s'f(-t)f'(-t) 

When  the  fitting  functions  are  restricted  to  trigonometric  functions  or 
polynomials,  0  approaches  zero  faster  than  f(-t)  grows,  since  B  <  1. 
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Thus,  F  will  reach  a  steady  state  value  after  a  sufficient  number  of  sam- 
pling intervals,  and  does  not  require  any  further  updating.   Once  steady 
state  has  been  reached  we  may  write 

F(t)  =  F(t-l)  =  FC-) 

If  one  of  the  fitting  functions  is  of  the  form  f(t)  =  e    ,  steady  state 
may  be  reached  only  if  B<e~  ^ .      In  the  case  of  trigonometric  and  polynomial 
models,  the  fastest  growing  function  will  be  the  highest  degree  polynomial 
term  present.   If  this  term  is  t  ,  then 

7  +  5.1n 


(1-3)°-" 


=  the  number  of  sampling  intervals  necessary  to  reach 
steady  state  (Brown,  2,  170), 
where  the  convergence  criterion  is  that  the  largest  quantity  added  to 
any  element  of  F  is  less  than  10   times  the  previous  value  of  that 
element . 

Having  reached  steady  state,  we  may  estimate  the  coefficients  by 

a(t)  =  F"^g(t) 

where  F  is  no  longer  time  dependent.   The  forecast  would  then  be  computed 
from 

i(t+T)  =  2'(t)f(T) 

=  [F"^g(t)]'f(T) 

=  g'(t)F"^f(T) 

=  g'(t)  c  (t) 
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where  c(t)  is  a  column  vector  dependent  on  the  values  of  the  fitting 
functions  at  time  i,  but  not  on  absolute  time.   Note  that  only  one  matrix 
inverse  will  be  necessary  as  long  as  B,  t,  and  the  fitting  functions  are 
not  changed,  so  the  computations  have  been  greatly  simplified.   When 
additional  data  is  available,  we  must  merely  update  g(t)  and  compute 
x(t+T). 

3.5  Derivation  of  General  Exponential  Smoothing 

We  have  said  that  F  will  converge  to  a  steady  state  and  may  be 
handled  as  a  program  constant  if:   (1)  successive  fitting  function  values 
may  be  generated  by  a  fixed  transition  matrix  L,  that  is,  f(t+l)  =  L  f(t); 
(2)  the  origin  of  time  for  the  fitting  functions  is  taken  as  the  most 
recent  observation.   Furthermore,  the  data  vector,  g(t),  may  be  defined 
recursively  from  the  previous  values  of  its  n  components,  g(t-l),  and  the 
new  observation.   We  wish  now  to  minimize  the  computational  effort  required 
in  obtaining  successive  estimates  of  the  coefficients,  a(t) . 

If  there  were  no  random  noise  in  the  data,  and  we  had  developed  a 
model  to  describe  the  process  exactly,  we  could  recursively  evaluate  a(t) 
from  a'(t-l).   In  this  case  x(t)  =  x(t),  so 

x(t)  =  I'(t-l)f(l)  =  x(t) 
x(t)  =  a'(t-l)f(l) 
=  a'(t-l)Lf(0) 
Note  also  that 

x(t)  =  a'(t)  f  (0) 

Equating  the  last  two  equations, 
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a'(t)£(o)  =  a'(t-l)L£(0) 
So. 

a^(t)  =  a^(l;-l)L  and  we  may  obtain  the  value  of  the  next 
observation,  that  at  time  t,  from 

i(t)  =  a'(t-l)  f(l) 
=  a'(t-l)L^(0) 
=  a'(t)f(0) 

We  cannot  expect  the  data  to  conform  exactly  to  the  model  of  the  process, 
however.   Each  observation  will  differ  from  the  corresponding  forecast 
by  an  amount  e , 

e(t)  =  x(t)  -  x(t). 

If  this  fitting  error  can  be  taken  into  account  in  successive  estimates 
of  the  coefficients,  we  may  still  obtain  a  recursive  relationship  for  a. 
In  order  to  account  for  these  fitting  errors,  we  shall  define  an  n  com- 
ponent smoothing  vector,  h,  which  depends  only  on  .f   and  the  fitting 
functions. 


h  = 


=  F  f (o)  =  the  smoothing  vector. 


h  may  also  be  treated  as  a  program  constant  for  a  particular  S  and  a 
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particular  model. 

As  in  the  previous  discussions,  the  minimum  discounted  squared 
residual  is  satisfied  when 

a(t)  =  F~^g(t) 

where  F~  is  used  in  its  steady  state.   This  may  also  be  written 

g(t)  =  Fa(t).  ■   ,• 

If  we  substitute  this  expression  into  the  recursive  relation  for  g(t) 
obtained  in  the  previous  section,  we  have 

Fa(t)  =  x(t)f(o)  +  eL'-'-Fa(t-l) 
or 

a(t)  =  x(t)F"-'-f(o)  +  6F"V-'-Fa(t-l) 

Note  that  the  coefficients  of  both  x(t)  and  a(t-l)  are  time  independent, 
and  that  F~  f (o)  is  h,  the  smoothing  vector.   By  manipulating  the 
equation  that  defines  F,  it  can  be  shown  that 

f'-'-l'-'t  =  L'  -  hf'(l) 

=  the  coefficient  of  a{t-l)  in  the  expression  for 
a(t)  above. 
Now  we  may  write 

a(t)  =  hx(t)  +  L'a(t-l)  -  hf '(l)a(t-l) 
But 
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f (l)a(t-l)  =  x(t), 
so, 

a(t)  =  hx(t)  +  L'a(t-l)  -  hx(t) 

-  L'a(t-l)  +  h[x(t)  -  x(t)] 
a(t)  =  L'a(t-l)  +  he(t). 

Thus  we  have  a  recursive  relationship  for  a  that  accounts  for  prtvious 
fitting  errors. 

3.6  Summary  of  General  Exponential  Smoothing 

The  relationships  used  are  valid  only  for  trigonometric  and  poly- 
nomial models,  and  exponential  models  satisfying  the  restriction  mentioned 
earlier  (Section  3.4).   The  general  forecast  equation  is 

x(t+T)  =  a'(t)f(T). 

The  fitting  function  coefficients  may  be  obtained  recursively  by 
a(t)  =  L'a(t-l)  +  he(t) 


where 


h  =  F'-'-f  (o). 


The  fitting  function  vector  nay  be  obtained  recursively  if  necessary  by 

f(t)  =  Lf(t-l). 

Appropriate  choices  of  3  and  the  initial  conditions,  a(o),are  important 
to  all  the  methods  discussed.   These  problems  will  be  discussed  in  the 
following  chapter. 
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4.   CHOOSING  THE  SMOOTHING  CONSTANT  AND  THE  INITIAL  COEFFICIENT  VECTOR 

A.l  Effect  of  Data  on  Discount  Factor  Choice 

There  are  essentially  three  ways  by  which  we  may  control  the  effici- 
ency with  which  our  forecasts  follow  the  actual  data: 

(1)  Choice  of  the  model 

(2)  Choice  of  the  discount  factor 

(3)  Choice  of  the  initial  conditions. 

The  first  control  was  discussed  at  length  in  an  earlier  chapter.   At  this 
time  we  wish  to  develop  some  guidelines  to  follow  in  choosing  6 ,  our  second 
control  over  the  model  efficiency.   The  literature  suggests  some  general 
rules  that  should  be  followed  in  making  a  B  choice,  however,  prior  research 
in  this  area  appears  to  be  inadequate.   Consequently,  much  of  the  research 
presented  in  this  thesis  is  based  on  the  effect  of  B  in  describing  the  true 
process. 

Initially,  our  choice  of  a  discount  factor  should  be  based  on  the 
nature  of  the  data  and  our  expectations  of  future  data.   If  the  deviations 
from  our  model  of  the  process  seem  to  be  mainly  random  fluctuations,  then 
a  large  B  value  will  probably  be  appropriate.   This  will  put  more  emphasis 
on  older  data  and,  hence,  tend  to  smooth  out  the  effect  of  the  fluctuations 
of  the  actual  data  on  our  forecasts.   If,  on  the  other  hand,  these  fluctu- 
ations are  large  in  magnitude  and  do  not  seem  to  imply  a  revision  of  the 
model,  a  small  3  will  put  greater  weight  on  more  recent  data,  and  allow  the 
model  to  react  more  quickly  to  these  fluctuations.   In  terms  of  the  fore- 
casting system  response  to  the  actual  data,  we  may  say: 
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(1)  For  large  values  of  S,  system  response  rate  is  low; 

(2)  For  small  values  of  6,  system  response  rate  is  high. 

If  the  investigator  expects  the  appearance  of  a  new  pattern  in  the  future 
values  of  the  actual  data,  i.e.,  if  he  expects  relatively  large  changes 
in  the  fitting  function  coefficient  values,  again  a  small  g  should  be  used 
to  allow  quick  response. 

We  shall  defer  discussion  of  actual  numerical  values  of  B  until  a 
later  section.   We  wish  first  to  develop  some  general  guidelines  on  the 
range  of  B,  as  dictated  by  the  data.   As  in  numerous  other  problems  we 
are  faced  with  two  conflicting  criteria  by  which  we  judge  the  efficiency 
of  our  S  choice,  that  is, 

(1)  Ability  to  respond  to  fluctuations  in  the  true  process 

(2)  Ability  to  smooth  out  random  fluctuations  in  the  data. 

A  smaller  6  gives  quicker  response  to  fluctuations  in  the  true  process, 
but  also  gives  more  drastic  response  to  the  random  fluctuations. 

4.2  Average  Age  of  the  Data 

Brown  (2,  116)  warns  that  if  you  must  use  a  B  smaller  than  0.7,  a 
larger  than  0.3,  you  may  actually  need  a  different  model,  specifically, 
a  constant  model.   He  feels  that  with  B  <  0.7  the  response  rate  is  so 
high  that  a  constant  model  may  be  quite  adequate.   It  may  be  that  an 
undetected  periodic  in  the  data  forces  the  use  of  a  high  S  .   Again  we 
may  use  a  more  reasonable  value  of  the  discount  factor  if  we  alter  the 
model  to  include  the  necessary  periodic  fitting  functions.   We  may  infer 
from  these  observations  that  it  would  be  wise  to  know  what  the  weighting 
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factors  will  be  for  specific  3  values.   Broun  (2,  107)  has  presented 
"average  age  of  the  data"  to  aid  the  investigator  in  grasping  the  meaning 
of  a  particular  13  value.   We  shall  apply  average  age  only  to  a  constant 
model  for  now.   In  Section  4.4  we  will  see  how  this  term  applies  to  more 
complicated  models.   We  may  define  the  age  of  the  data  as  being  one  for 
the  present  observation,  two  for  the  previous  observation,  three  for  the 
one  before  that,  and  so  on.   So,  the  observation  N  sampling  intervals  age 
will  have  age  N.   If  we  multiply  the  age  of  each  data  value  used  in  fore- 
casting by  the  weight  received  by  that  data,  we  have  a  figure  representing 
the  average  age  of  all  the  data  used.   Hence,  for  a  moving  average  of 
period  N, 

k  =  i  (0)  +  i  (1)  +  i  (2)  +  .  .  .  +  i  (N-1) 

k  =^^ 
2 

where 

k  =  average  age  of  the  N  data  points 

—  =  the  weight  received  by  each  data  point 

Applying  this  to  exponential  smoothing,  we  see  that  the  weight  given  to 

data  n  sampling  intervals  old  is  ag  .   Hence,  the  average  age  of  the 

data  would  be 

k  =  a(o)  +  aB(l)  +  aB^(2)  +  .  .  .  +  a6^(T) 

T 

=  a   I  Bn  " 

n=0 
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where 

T  =  the  total  number  of  observations  used  in  making  the 
forecast. 

If  T  is  large,  that  is,  if  we  have  much  data  from  which  to  make  a  fore- 

T  ~ 

cast,  B  will  be  very  small.   We  may  then  approximate  k  by  letting  T 

approach  infinity  without  introducing  significant  error  into  the  calcu- 
lations (Brown,  2,  135).   So, 


a  E  n^ 
n=0 


'-^ 


"•(1-3) 

a6 
a 

k  =  B/a  . 

We  might  want  to  define  an  exponential  smoothing  process  as  being  equiv- 
alent to  a  moving  average  with  the  same  average  age  of  the  data.   That 
is. 

S   N  -  1 


N  -  1 


N  +  1 

1  +  B   1  + 
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This,  of  course,  would  only  apply  to  exponential  smoothing  with  a  constant 
model  of  the  process. 


4.3  System  Response  to  Standard  Input  Signals 

We  may  gain  some  feeling  for  the  effect  of  3  on  the  system  by  analyzing 
the  system  response  to  various  standard  input  signals.   The  inferences  to 
be  dratjn  are  only  of  a  general  nature,  since  the  constant  process  model 
assumed  here  would  often  be  inadequate. 

Before  we  study  system  response,  it  is  necessary  to  establish  a  few 
definitions.   Any  discrete,  linear,  time-invariant  system  can  be  completely 
described  by  its  impulse  response  (Brown,  2,  112).   All  of  our  discussions 
have  been  based  on  the  assumption  of  discrete  data,  that  is,  data  taken 
only  at  specific  and  equal  time  intervals.   Linear  implies  the  following 
system  characteristic  regarding  input  to  and  output  from  the  system; 
If  input  X- (t)  produces  output  y. (t) ,  and  input  x„(t)  pro- 
duces output  y^(t),  then  input  ax  (t)  +  bx„(t)  will  produce 
output  ay  (t)  +  by„(t),  where  a  and  b  are  constants. 
The  system  is  time-invariant  if  a  particular  input  signal  will  cause  the 
same  system  output  at  any  time,  t. 

A  unit  impulse  signal  may  be  as  shown  in  Figure  4:3:1. 


1.0 


x(t) 


Fig.  4:3:1  Unit  impulse  function 
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It  may  be  written  as 

x(t)  =  1         t  =  0 


x(t)  =  0         t  j«  0 


where 


x(t)  =  the  input  signal 

t  =  time. 
Since  we  have  not  defined  our  exponential  smoothing  system  for  negative 
time,  we  may  assume  a  time  range  of  0  <t  <  ".   Let  h(t)  be  the  unit  impulse 
response  of  the  system  t  periods  after  the  unit  impulse  input.   The 
response  of  a  linear  system,  y(t),  to  an  arbitrary  signal  x(t)  may  be 
expressed  as 

y(t)  =   I  x(n)h(t-n) 

n=-«> 

Since  our  system  is  undefined  for  t<0. 


y(t)  =  Z     x(n)h(t-n) 
n=0 


=  x(t)*h(t) 

where  the  symbol  *  implies  a  convolution  of  the  sequence  [x(t)]  with  the 
sequence  [h(t)].   Note  from  the  previous  equation  that  for  a  linear 
function  of  past  data,  as  in  exponential  smoothing  with  a  constant  model, 


h(t-n)  =  aw(t-n)  =  a3  . 


Hence, 
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y(t;)  =  a   I     x(n)w(t:-a)  . 

k=0 

So,  having  found  an  expression  for  the  unit  impulse  response  for  exponen- 
tial smoothing,  ue  may  find  the  system  response  to  an  arbitrary  input  x(t) . 

First,  let  us  simplify  the  computations  by  transforming  the  convolu- 
tion equation  above.   The  Z-transform  (Appendix  D)  is  particularly  conven- 
ient for  this  analysis,  since  the  convolution  becomes  a  multiplication  in 
the  Z-doraain.   The  Z  -  transform,  F(Z),  may  be  written 

F(Z)  =  I    f(n)z'' 

n=-«> 

where 

f(n)  =  the  function  to  be  transformed 
F(Z)  =  the  Z  -  transform  of  f(n). 

Note  that  the  previous  equation  for  y(t)  may  now  be  expressed  as  a  multip- 
lication, 

Y(Z)  =  X(Z)'H(Z) 
where 

Y(Z)  =  the  Z-  transform  of  y(t) 

XCZ)  =  the  Z-  transform  of  x(t) 

H(Z)  =  the  Z-  transform  of  h(t). 

Once  we  have  found  Y(Z),  we  may  apply  an  inverse  transform  to  obtain  y(t) . 
The  system  response,  y(t) ,  is  determined  below  for  several  standard  inputs. 

For  a  unit  impulse  input  at  time  zero,  as  described  above,  we  may 
write 
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Now, 


So, 


X(2)  =  E  x(n)Z 

n=0 


l'Z°  +  0  +  0  +  .  . 


-  1 
H(Z)  =  Z     aB"z"  =  a  E   e"z" 


n=0 


"^l4^^ 


n=0 


Y(Z)  =  X(Z)-H(Z) 


^  •  '='^1::^^ 


Y(Z)  =  !•  !:.  aeV 

n=u 


y(t)  =  aS>\ 


Hence,  y(t)  for  a  unit  impulse  input  is  precisely  the  weight  given  to 
the  data.   Note  that  the  system  response  is  dependent  upon  our  choice  of 
£,  and  hence  of  a.   This  effect  may  be  seen  in  Figure  4:3:2. 


x(t),y(t) 


1  n 

■ x(t) 

) 

\~ 

— y(t) 

n 

I 

1          1 

^^^~T--^ 

-2    -1     0     1     2     3     A     5 
Fig.  4:3:2  System  response  for  unit  impulse  input  to  a  constant  process 
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For  lower  values  of  3  the  system  will  adjust  back  to  its  original  state 
at  a  faster  rate,  since  w(t-n)  =  3  will  decrease  at  a  faster  rate  than 
it  would  for  a  higher  6  value. 

Suppose  the  input  signal  is  a  unit  step,  as  shown  in  Figure  4:3:3, 


1.0 


x(t) 


J I I I L. 


Fig.  4:3:3  Unit  step  function 

This  may  be  written 

x(t)  =0        t  <  0 
x(t)  =1        t  >  0. 

Determining  the  Z  -transforms 

x(z)=X  ^"•^"^ilz 


H(2) 


1-6Z  ' 


We  now  have 


Kow, 


^(^^  =  (i=z)-^iril) 


I  y(n)z"  =  C  !:   z")-(  Z  agV) 
n=0  n=0      n=0 


I  aZ   £  B 


n=0     k=0 
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and 


y(t)  =  1 


^  tt(l-e    )  ,n 
n=0  "  ^-3 


t+1 


Note  again  the  effect  of  B  on  y(t).   For  lower  S  values  the  system  will 
adjust  to  the  input  at  a  much  faster  rate.   This  effect  is  shown  in 
Figure  4:3:4. 


30. 


Periods   20 
to  rise 
to  90% 
of  step.  10 


1.0 


Fig.  4:3:4  Effect  of  discount  factor  on  system  response  to  a  unit 
step  input 

Since  exponential  smoothing  gives  some  weight  to  all  past  data,  it  will 

never  actually  reach  the  step,  but  approaches  it  asymptotically  as 

shown  in  Figure  4:3:5. 
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1.0. 


K(t),  y(t)' 


y(t:) 


-. x(t) 


Fig.  4:3:5  System  response  to  unit  step  input 

Through  the  same  procedure  we  will  find  that  for  a  ramp  input  signal 

x(t)  =  t»     t_>  t, 
y(t)  =  t  -^  d-s'^). 

Again  we  see  the  dependence  of  system  response  on  our  choice  of  the  dis- 
count factor. 

Summarizing  our  observations  of  system  response  characteristics, 
let  us  re-iterate  once  more  the  general  rule  regarding  our  discount 
factor  choice: 

The  lower  our  choice  of  B,  the  higher  will  be  the  system  response  rate. 


4.4  Equivalent  Discount  Factor 

Ue  may  describe  an  exponential  smoothing  process  as  having  an  order 
equal  to  the  number  of  coefficients  in  the  corresponding  model  or  the 
number  of  coefficients  to  be  estimated.   The  constant  process  model  with 
which  we  have  been  dealing  represents  first  order  exponential  smoothing. 
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Let  us  make  the  safe  assumption  that  smoothing  of  any  order  n  includes 

,    th    , 
a  constant  process  term  in  the  model.   There  is  some  g   tor  n   order 

smoothing  that  will  give  the  same  estimate  for  the  constant  term  as  S^ 

in  single  smoothing  with  a  constant  model. 

Consider  the  response  of  a  system  to  a  unit  impulse  input  at  time 

zero.   The  estimate  of  the  constant  term  coefficient,  a  (t) ,  is  given 

by 

;,(t)  =  1  -  b;; 

for  exponential  smoothing  of  order  n  (Brown,  2,  147).   We  may  now  equate 
the  a  (t)  estimates  for  single  smoothing  and  smoothing  of  order  n. 

;^(t)  =  1  -  B^  =  1  -  b;; 

Hence,  the  equivalent  discount  factor  is  defined  by 

B^  =  b"  . 
1    n 

Thus,  if  B  is  the  discount  factor  used  in  an  n   order  exponential 
n 

smoothing  process,  a  S  equal  to  B  would  give  the  same  estimate  of  a^(t) 
in  single  smoothing.   Stated  another  way,  if  S.  is  the  discount  factor 
used  in  single  smoothing,  a  B  equal  to 


B  =Bt^" 
n    1 


will  give  the  same  a  (t)  in  n   order  exponential  smoothing. 
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We  infer  from  this  discussion  Chat  an  n   order  snoothing  system 
will  have  response  characteristics  similar  to  a  first  order  system  with 
an  equivalent  discount  factor.   Let  us  reiterate  the  major  inference  to 
be  drawn  here: 

A  smaller  B  gives  more  weight  to  more  recent  data  than  does  a  larger 

$  value, 

4.5  Choosing  the  Initial  Coefficient  Vector 

Exponential  smoothing  always  requires  some  initial  condition  from 
which  to  start  the  smoothing  process.   This  logically  entails  some  initial 
estimate  of  the  coefficients  of  the  model  fitting  functions,  or  a(0). 
Again  we  must  note  that  very  little  conclusive  research  on  the  choice  of 
a(0)  was  found  in  the  literature.   The  relative  merit  of  an  accurate  choice 
is  a  subject  for  debate,  however,  since  the  need  for  accuracy  here  is 
greatly  dependent  upon  our  choice  of  the  discount  factor  and  the  amount  of 
past  data  available.   For  example,  if  we  have  very  little  past  data  at  the 
time  the  smoothing  process  begins,  we  may  want  to  estimate  a(0)  as  accurately 
as  possible  to  shorten  the  length  of  time  necessary  for  the  smoothing 
system  to  stabilize  in  the  light  of  new  data.   On  the  other  hand,  if  we 
have  ample  past  data,  the  smoothing  system  will  have  more  information  from 
which  it  may  stabilize  before  making  a  forecast.   The  dependence  of  a(0)  on 
3  may  be  expressed  generally  in  this  v;ay: 

If  you  are  very  confident  that  your  estimate  of  a(0)  is  a  good 
one,  use  a  large  3,  since  quick  system  response  to  errors  in  a(0) 
will  probably  not  be  necessary.   Conversely,  if  you  have  very 
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little  confidence  in  your  initial  estimate,  use  a  smaller  S, 

so  the  system  may  adjust  to  the  possible  error  more  quickly. 

There  are  several  means  of  obtaining  the  initial  coefficient  estimate. 
If  the  investigator  is  thoroughly  familiar  with  the  process  Involved,  he 
may  be  satisfied  with  a  subjective  estimate  of  a(0).   He  may  employ  a 
graph  of  past  data  against  time  to  estimate  the  initial  conditions.   Perhaps 
the  most  logical  approach  is  to  use  the  least  squares  regression  coefficients. 
This,  however,  requires  some  past  data,  the  more  the  better.   If  we  have 
no  past  data,  we  may  be  able  to  compare  the  process  with  another  one  for 
which  we  have  ample  data.   This  might  be  true  in  a  case  where  x(t)  repre- 
sents demand  for  some  new  product.   Past  data  on  a  similar  product  may 
serve  as  a  good  predictor  of  a(0)  for  the  new  product. 

Brown  (2,  194)  presents  a  scheme  for  estimating  the  initial  coeffic- 
ient values  for  models  containing  a  constant  and  some  periodic  terms. 
First,  if  the  constant  coefficient  is  a  (t) ,  let  a  (0)  be  the  average  of 
the  past  data  over  any  number  of  full  cycles.  Let  a„(t)  and  a_(t)  be 
the  coefficients  of  a  pair  of  periodic  terms  for  a  period  of  p  defined  by 
the  fitting  function  pair  for  which  we  want  coefficient  estimates.   If 
there  are  several  cycles  of  past  data  available,  average  corresponding 
data  samples  across  them.   If  there  are  no  data,  estimate  the  cyclic 
pattern.   By  one  of  these  methods  we  obtain  a  sequence  of  p  numbers 

X, ,  x_,  .  .  .,  X  .   The  initial  coefficient  estimates  are  then  obtained 
12         p 

from 
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where 


and 


a„(0)  =  —  i;  X  sin  ^^  k. 
2      p  k=l  k    p 


2  P  -    2Tr 
3      P  k=l  '^    P 


p  =  the  number  of  sampling  intervals  in  the  period  under 
study 

k  ranges  from  1  to  p. 
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3.   GATHERING  THE  DATA 

5.1  Time  Series  Requirements 

The  general  exponential  smoothing  techniques  developed  earlier 
require  the  formulation  of  some  model  that  is  representative  of  the  tlirie 
series  for  which  we  want  a  forecast.   Generally,  the  forecaster  cannot  ade- 
quately determine  the  fitting  functions  to  be  included  in  this  model  with- 
out some  analysis  of  past  data.   Data  analysis  may  also  aid  the  forecaster 
in  making  a  good  initial  coefficient  vector  choice.   Such  analyses  were 
treated  in  previous  chapters. 

The  past  time  series  data  must  conform  to  the  analyses  for  which  it 
will  be  used.   So,  it  should  be  discrete  data,  taken  at  specific  and 
equal  sampling  intervals,  from  the  present  time  back  to  some  past  time. 
Furthermore,  the  data  should  be  in  a  form,  i.e.,  dollars,  gallons,  etc., 
that  is  easily  manipulated  and  Interpreted.   Consistency  of  the  data  is 
vital  in  that  a  particular  time  series  value  should  truly  represent  the 
time  series  process  for  some  particular  sampling  interval. 

Data  exhibiting  such  characteristics  is  not  always  readily  available. 
Ideally,  a  forecaster  would  want  to  design  the  record-keeping  system 
from  which  he  will  obtain  the  time  series.   Unfortunately,  this  is  not 
always  possible,  and  we  must  use  what  data  is  available  at  the  present 
time.   There  are,  however,  some  alternatives  to  the  manner  in  which  the 
data  Is  manipulated.   It  may  be  desirable  to  transform  the  data.   By  doing 
so  we  may  produce  a  more  meaningful  series,  or  one  that  is  more  easily 
predicted.   Transformation  of  data  can  usually  be  performed  by  some 
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mathematical  operation,  such  as  a  multiplier,  log,  etc. 

Suppose  the  time  series  for  which  we  want  forecasts  is  incomplete, 
or  not  always  available.   We  may  find  that  there  is  some  related  time 
series  that  may  be  more  easily  predicted.   If  there  is  some  functional 
relationship  between  the  two  series,  we  can  obtain  the  desired  forecast 
from  a  forecast  of  the  related  series. 

Another  possibility  is  sub-classification  of  the  time  series.   It 
may  be  easier  to  forecast  separate  parts  of  a  time  series  and  sum  these 
forecasts  to  get  a  forecast  of  the  whole.   This,  of  course,  would  involve 
more  work  than  a  single  forecast,  but  the  need  for  accuracy  may  warrant 
such  a  procedure. 

In  our  discussion  of  model  types  (Section2.$  ,  we  implied  a  further 
forecasting  alternative.   We  may  be  able  to  find  one  or  a  combination  of 
other  series  that  acts  as  a  leading  series  to  the  one  which  we  want  to 
predict.   That  is,  this  other  series  may  serve  as  an  indicator  of  what  the 
original  series  will  look  like  a  certain  number  of  sampling  intervals  in 
the  future.   In  this  case  we  may  wish  to  use  an  autoregressive  model  in 
the  adaptive  smoothing  technique  mentioned  earlier. 

5.2  Actual  Data  Used 

The  data  used  in  this  research  consists  of  four  time  series  taken 
from  the  Manhattan  Wholesale  and  Retail  Floral  Companies  of  Manhattan, 
Kansas.   These  time  series  represent  demand  for  four  types  of  plants 
produced  by  Manhattan  Wholesale: 

1.   Potted  chrysanthemums  sold  only  through  Manhattan  Ifljolesale. 
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2.  Potted  chrysanthemums  sold  only  through  Manhattan  Retail. 

3.  Potted  poinsettias  sold  through  Manhattan  Wholesale. 

4.  Potted  lilies  sold  through  Manhattan  Wholesale. 
Chrysanthemums  are  produced  and  sold  throughout  the  year.   The  two  time 
series  related  to  this  product  group  represent  end-of -month  dollar  receipts 
on  sales  of  these  plants.   The  latter  two  products  are  sold  strictly  on 

a  seasonal  basis.   Poinsettias  are  sold  only  as  a  Christmas  season  item, 
while  lilies  are  sold  only  as  an  Easter  plant.  Data  on  seasonal  items  is 
in  number  of  pots  sold  per  day. 

Both  the  chrysanthemum  time  series  extend  over  a  period  from  September, 
1958,  to  May,  1957,  inclusive.   This  gives  a  total  of  105  observations,  or 
more  than  eight  complete  yearly  cycles.   Figures  reflect  the  actual 
receipt  of  payments  during  a  month  for  chrysanthemum  sales.   So,  a  figure 
does  not  necessarily  represent  the  sales  that  occurred  during  a  particular 
month,  since  not  all  buyers  pay  their  bills  promptly  each  month.   Neverthe- 
less, this  was  the  only  data  available  on  these  sales  at  the  time  of  this 
research.   This  fact  may  explain  the  high  noise  level  of  the  data,  as  we 
shall  see  later.   Plots  of  the  data  are  shown  in  Appendix  E. 

It  should  be  noted  here  that  the  data  available  at  the  time  a  fore- 
cast is  desired  will  probably  be  anything  but  ideal.   Hence,  future  data 
compilation  should  be  altered,  if  necessary,  to  best  reflect  the  character- 
istics of  the  data  for  which  a  forecast  is  desired. 

The  retail  receipts  chrysanthemum  time  series  shows  some  growth. 
This  can  probably  be  attributed  to  a  combination  of  inflation  and  market 
expansion.   However,  the  market  in  the  case  of  retail  sales  is  probably 
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nearly  saturated,  since  it  is  confined  mainly  to  Manhattan  alone.   The 
wholesale  chrysanthemum  time  series  shows  a  much  greater  growth  trend. 
Most  of  this  is  due  to  market  expansion.   The  wholesale  market  at  present 
includes  numerous  cities  other  than  Manhattan.   Though  many  of  the  cities 
to  which  deliveries  now  go  are  near  saturation,  the  market  as  a  whole  is 
not. 

The  poinsettia  time  series,  also  plotted  in  Appendix  E,  covers  the 
demand  periods  for  the  years  1963  through  1966.   Thus,  four  complete  sets 
of  data  were  available.   The  market  area  is  roughly  the  same  as  for 
wholesale  chrysanthemum  sales.   The  demand  period  per  year  for  poinsettias 
was  arbitrarily  set  as  the  thirty-five  days  immediately  preceding  Christmas 
Day.   Since  the  earliest  purchase  date  for  the  data  used  was  thirty-three 
days  before  Christmas,  all  data  is  included  within  the  range  of  the  thirty- 
five  day  period,  and  we  have  a  total  of  140  observations.   Normally, 
deliveries  were  not  made  on  Sundays,  so  the  data  shows  a  fluctuation  about 
the  zero  line,  unlike  the  two  series  previously  discussed. 

The  lily  time  series.  Appendix  E,  has  a  form  similar  to  that  for  the 
poinsettias,  as  might  be  expected.   The  data  also  covers  the  sales  periods 
for  the  years  1963  through  1967.   This  time,  hoxv'ever ,  a  demand  period  of 
fourteen  days  before  Easter  is  large  enough  to  include  all  of  the  data. 
Hence,  we  have  70  observations. 

The  market  growth  may  be  seen  in  these  seasonal  items,  poinsettias 
and  lilies,  in  the  form  of  an  amplitude  growth  rather  than  a  trend  line. 
For  these  items  the  data  available  was  more  complete;  that  is,  it  was 
broken  down  by  customer.   Hence,  we  are  able  to  segregate  the  market  for 
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poinsettlas  and  lilies  in  some  fashion,  if  desirable.   By  comparing  the 
total  number  of  potential  buyers  (usually  retail  florists)  in  each  city  with 
the  number  to  whom  we  now  deliver,  the  local  market  saturation  becomes 
apparent.   We  might  interpret  this  to  mean  that  unless  new  market  areas  are 
exploited,  the  growth  of  sales  may  level  off  somewhat. 

Since  all  of  the  customer  outlets  are  in  essentially  the  same  geo- 
graphic area,  segregation  by  geographic  area  was  not  feasible.   Though 
most  cities  in  the  market  area  are  saturated,  there  are  one  or  two  larger 
cities  in  which  relatively  few  outlets  buy  from  Manhattan  Floral  l-Jhole- 
sale  Company.   Hence,  segregation  of  customers  by  city  size  is  meaning- 
less.  Data  was  not  available  by  which  segregation  by  customer  size  could 
be  performed. 

Another  point  should  be  clarified  before  we  continue.   The  existing 
productive  capacity  of  the  Manhattan  Floral  Wholesale  Company  is  nearing 
full  utilization,  so  much  further  market  expansion  would  be  impossible 
without  an  expansion  of  productive  capacity  also. 

The  possibility  of  correlating  some  independent  variable  or  variables 
with  a  particular  time  series  was  also  looked  at  briefly.   The  obvious 
problem  here  is  finding  a  time  series  that  has  a  high  correlation  with 
one  of  the  afore  mentioned  series,  and  is  also  known  for  future  time. 
In  addition,  assuming  a  related  series  is  available  we  incur  the  added 
expense  of  keeping  multiple  data  records.   Such  a  series,  though  none 
could  be  found,  would  have  enabled  the  application  of  adaptive  smoothing 
with  an  auto-regressive  model. 

Due  to  the  data  limitations  presented  above,  all  time  series  were 
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used  in  the  form  in  shich  they  were  gathered.   The  following  chapter  con- 
tains the  results  of  the  application  of  general  exponential  smoothing  to 
these  four  time  series. 
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6.   DISCUSSION  OF  RESULTS 

6.1  Experimental  Procedure 

All  four  time  series  were  initially  plotted  (Appendix  E) .   The  forms 
of  the  resulting  curves  were  noted  in  order  to  gain  some  insight  into  the 
fitting  functions  that  should  be  included  in  the  models  for  each  time 
series  process.   In  addition,  detrending,  auto-correlation>  and  harmonic 
analyses  were  performed  on  each  time  series.   From  the  results  of  these 
analyses  and  from  intuition,  various  models  were  proposed.   These  models 
were  then  used  in  the  general  exponential  smoothing  techniques  for  various 
initial  conditions,  various  initial  fitting  function  coefficient  vectors, 
and  a  range  of  discount  factor  values. 

In  applying  general  exponential  smoothing  to  a  time  series,  the 
steady-state  F-matrix  was  first  determined,  for  the  particular  ^  and  a(o), 
from  the  data  given.   Then,  using  the  methods  discussed  in  Section  3.5, 
forecasts  were  made  for  all  sampling  intervals  for  which  data  was  actually 
available.   So,  starting  at  time  zero,  a  forecast  was  made  for  the  first 
time  series  value  and  compared  to  the  actual  value  at  observation  one. 
Based  on  this  comparison,  a  forecast  was  then  made  for  time  series  value 
two,  and  compared  to  actual  value  two,  and  so  on.   Each  forecast  was  pro- 
jected one  sampling  interval  into  the  future.   Forecasts  were  performed 
only  for  time  series  values  for  which  actual  data  existed.   In  this  way 
we  could  measure  the  relative  efficiency  of  a  particular  model  -  initial 
condition  -  discount  factor  combination.   Note  that  it  would  be  impossible 
to  measure  the  accuracy  of  a  forecast  for  which  no  actual  data  existed. 
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A  discussion  of  the  various  efficiency  criteria  follows. 

Generally,  the  experimental  procedure  for  a  single  time  series  may 
be  summarized  as  follows: 

(1)  Plot  the  time  series  against  time. 

(2)  Perform  a  detrending  analysis. 

(3)  Perform  an  auto-correlation  analysis. 

(4)  Perform  a  harmonic  analysis. 

(3)   Formulate  various  models  to  be  tried. 

(6)  Determine  the  F-matrix, 

(7)  Use  general  exponential  smoothing  to  forecast  the  original 
time  series.   Do  this  for  various  values  of  a(o)  and  S. 

(8)  Compute  the  efficiency  criteria  for  each  exponential  smoothing 
iteration. 

(9)  Compare  the  results  and  draw  conclusions. 

6.2  Efficiency  Criteria 

In  order  to  evaluate  and  compare  the  forecasting  efficiency  of  the 
models  used,  it  is  necessary  to  establish  certain  efficiency  measurement 
methods.   The  obvious  approach  to  this  problem  is  to  measure  the  accuracy 
of  fit  of  the  general  exponential  smoothing  process  using  a  particular 
model- initial  condition-discount  factor  combination.  With  this  goal  in 
mind,  the  following  efficiency  criteria  were  proposed: 

9 

(1)  Variance  of  the  forecasts,  a 

(2)  Standard  deviation  of  the  forecasts,  a^ . 

(3)  Mean  error,   e 


(4)  Mean  absolute  error, |e| 

(5)  Mean  percent  error,   e%. 

—2 

(6)  Variance  of  the  error,  oc  . 

(7)  Error  standard  deviation  about  the  mean  error,  o— . 

2 
The  first  criterion  above,  0_,  is  not  truly  the  variance  of  the 

forecasts,  but  the  variance  of  the  fitting  errors  about  zero.   The  compu- 
tations necessary  to  determine  the  true  forecast  variance  are  quite 
lengthy,  so  this  simple  computation  has  been  used  instead  (Brown,  2,  393). 
We  shall  call  the  former  quantity  a  forecast  variance  to  distinguish  it 
from  the  error  variance  about  the  mean  error,  which  is  defined  later. 
So,  let  us  compute  the  variance  of  the  forecasts  as: 

N 

2         2 

0^  =  I      el  /N 

1=1 
where 


X.  =  the  i   actual  observation 

X.  =  the  exponential  smoothing  forecast  of  the  i 

observation 
N  =  the  total  number  of  data  points  observed  and  forecasted. 

Similarly,  we  may  now  estimate  the  standard  deviation  of  the  forecasts  as: 


^2- 


The  mean  error  may  be  computed  from 
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N 
=  E  e./N. 
i=l  ^ 


Note  that  e  may  be  positive  or  negative.   If  we  consider  a  shortage  and 
a  surplus  as  being  equally  detrimental,  we  may  prefer  to  look  at  the 
mean  absolute  error, jel,  or 

171=  I   1..I/N. 
i=l 

A  similar  error  measurement  method  is  the  mean  per  cent  error,   e%,  or 


e%  =  100-  I  ^^^  I   N. 
i=l   i 


The  difference  between  this  efficiency  criterion  and  the  mean  absolute 
error  needs  emphasis.  |el  is  in  the  same  units  as  the  original  time  series, 
whereas  e  %  is  a  percentage  of  the  actual  data.   Thus,  when  |  £  .1  is  small, 
£,%  may  be  quite  large  if  x.  is  also  small.   The  converse  is  also  true. 
The  variance  of  the  error  is  a  true  variance  computed  relative  to 
the  mean  of  the  errors.   It  may  be  written: 


2  -2 

a-  =  I   (e.  -  e)  /N 

^         i=l   "■ 

Then  a—  is  merely  the  square  root  of  the  variance  above,  lio— .   Note  that 
if  there  were  no  noise  in  the  data,  and  if  our  model  truly  described  the 
process,  then  0"— would  be  identical  to  a  ,  since  e  would  be  zero. 
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Not  all  of  the  criteria  above  were  studied.   Since  variance  and 
standard  deviation  basically  represent  the  same  error  characteristic, 
only  standard  deviation  was  analyzed.   Thus,  we  are  left  with  five 
efficiency  criteria  with  which  to  compare  various  models: 

(1)  Standard  deviation  of  the  forecasts. 

(2)  Mean  error, 

(3)  Mean  absolute  error. 
(A)  Mean  percent  error. 

(5)   Error  standard  deviation  about  mean  error. 
Mean  error  proved  to  be  difficult  to  interpret  and,  consequently,  it  was 
discarded  as  an  efficiency  criterion.   The  mean  error  over  the  entire 
range  of  available  data  represented  such  a  small  percentage  of  the  corres- 
ponding time  series  average  that  a     and  0—  were  nearly  identical.   Hence, 
the  standard  deviation  of  the  forecasts  was  also  discarded  as  an  efficiency 
criterion. 

6.3  Computer  Programs 

All  computations  were  performed  on  the  IBM  System  360/50  digital 
computer  at  Kansas  State  University,   The  programs  used  are  shown  in 
Appendix  F.  The  most  important  of  these  programs  are: 

(1)  AUTODET 

(2)  POWSPEC 

(3)  PLOTOOl 

(4)  EXPOSMO 

AUTODET  first  detrends  the  original  time  series  using  a  first  degree 
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polynomial.  It  then  computes  the  auto-correlation  coefficients  of  the 
detrended  data  for  all  lags  specified  by  the  user.  POWSPEC  uses  the 
detrended  data  from  AUIODET  to  perform  a  harmonic  analysis  on  the  time 
series.  It  computes  values  of  a  periodogram  for  the  periods  specified 
by  the  user.  The  autocorrelation  function  and  periodogram  may  then  be 
plotted  using  PLOTOOl. 

EXPOSMO  is  the  computer  program  which  performs  the  general  exponen- 
tial smoothing  process.   Inputs  to  the  program  include  the  original 
time  series,  the  initial  coefficient  vector,  the  transition  matrix,  and 
the  values  of  the  discount  factor  the  user  wishes  to  try.   The  program 
computes  and  lists  the  transition  matrix,  the  F  and  F- inverse  matrices, 
the  smoothing  vector,  the  variances  of  the  coefficients  in  the  model, 
the  forecasts,  and  the  error  criteria  values.   In  addition,  the  program 
outputs  a  plot  against  time  of  the  actual,  forecast,  and  error  data. 

EXPOSMO,  AUTODET,  and  POWSPEC  are  taken  from  programs  written  by 
Mr.  Raymond  Miller,  a  former  student  at  Kansas  State  University.   For 
a  more  complete  explanation  of  their  operation  refer  to  his  thesis (Miller , 
7,  65).   Only  minor  changes  were  made  to  AUTODET  and  POWSPEC.   Significant 
changes  were  made  to  the  input-output  structure  of  the  exponential 
smoothing  program  written  by  Miller  in  developing  EXPOSMO.   The  basic 
computational  procedure,  however,  was  not  altered. 

6.4  Effect  of  the  Initial  Coefficient  Vector 

Two  approaches  to  choosing  the  initial  coefficient  vector  were  tried. 
In  one  approach  the  model  coefficients  determined  in  the  least  squares 
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detrending  process  were  used.   Recall  that  in  detrending  we  used  first 
degree  polynomial  regression,  so  coefficients  were  computed  only  for  the 
two  corresponding  polynomial  fitting  functions.   In  the  initial  coeffic- 
ients were  used  when  applicable;  all  other  initial  coefficients  were  set 
at  zero. 

The  other  choice  of  initial  conditions  was  purposely  made  to  be 
erroneous;  that  is,  initial  coefficients  were  purposely  set  far  from 
their  proper  values.   The  resulting  exponential  smoothing  fitting  errors 
were  very  large  for  earlier  sampling  intervals. 

Efficiencies  around  the  best  value  of  the  discount  factor  (the  6 
which  minimizes  the  efficiency  criteria  for  that  model)  are  similar  for 
either  initial  coefficient  vector  choice;  however,  for  the  second  choice 
the  degree  of  error  increases  very  rapidly  as  we  move  away  from  this 
best  B.   Thus,  when  a  bad  choice  of  initial  conditions  is  made,  the 
system  is  less  stable,  making  the  proper  choice  of  the  discount  factor 
more  critical. 

6.5  Model  Identification  Notation 

Throughout  the  remainder  of  this  text,  it  will  be  convenient  for  us 
to  develop  a  shorthand  notation  for  identification  of  the  various  models 
used.  Basically,  for  all  models  the  fitting  functions  used  fall  into  one 
of  three  general  categories: 

(1)  Cyclic  terms 

(2)  Polynomial  terms 

(3)  Combination  of  cyclic  and  polynomial  terms. 
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Cyclic  terms  are  sine  and  cosine  function  pairs  for  some  specific  integral 
period.   We  shall  denote  such  terms  as  Cp,  where  C  tells  us  that  both  a 
sine  and  cosine  fitting  function  are  included  and  p  is  the  number  of 
sampling  intervals  taken  as  the  period  for  these  functions.   Thus,  the 
notation  C 12  would  indicate  that  the  following  terms  are  included  in  the 
model ; 

^.  .ZTTt  ,     „   2Trt 

If  more  than  one  pair  of  cyclic  terms  is  present,  they  will  be  written 
Si'  ^2'  Pn'  °' 

a.Sin^  +  a.^  cos^  +  a.^,Sir^ 
1    pi    1+1    pi    i+2    p2 

+  a  .  ,  „Cos — r-  +  a  ,  , ,  Sinr—  +  a  .  ,  ,Cos . 

1+3    p2    1+4    pn    1+5    pn 

Polynomial  fitting  functions  include  all  terms  of  the  form  shown 
in  Section  3.3.   In  all  of  the  models  containing  polynomial  terms,  a 
complete  polynomial  of  some  degree  is  used.   Note  that  a  polynomial  of 
degree  zero  designates  the  inclusion  of  a  constant  fitting  function  only. 
We  shall  denote  polynomials  as  Pd ,  where  P  stands  for  polynomial  and  d 
is  an  integer  designating  the  degree  of  the  polynomial  used.   Thus,  the 
notation  P3  tells  us  that  the  following  terms  are  present  in  the  model: 

a^  +  a^^t  +^  t(t-l)  +f2  t(t-l)(t-2). 
6 

By  combinations  of  cyclic  and  polynomial  terms,  we  mean  the  product 
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of  such  fitting  functions.   These  combination  fitting  functions  are 
encountered  in  only  one  form  In  this  research,  that  Is,  the  product  of 
a  cyclic  pair  and  time,  t.   Since  such  terms  produce  linearly  increasing 
values  for  the  amplitude  of  the  sine  and  cosine  functions,  we  may  refer 
to  them  as  amplitude  growth  terras,  or  more  simply,  growth  terms.  We 
shall  designate  such  terms  with  the  notation  Gp,  where  G  indicates  growth 
terms  and  p  is  the  period  of  the  cyclic  terras  to  which  this  linear  ampli- 
tude growth  is  applied.   Thus,  G6  tells  us  that  the  following  terms  are 
included  in  the  process  model: 

>o-   2Tit  ,       ^„       2TTt 

a  .  tSm  —7—  +  a  .  , ,  tCos— 7—  . 
1      6     1+1     6 

Models  containing  terms  from  more  than  one  of  the  categories  discussed 
above  may  be  designated  in  the  same  fashion,  with  commas  separating  the 
fitting  function  notations  of  different  categories.   So,  the  notation  C6, 
12,  P2,  G12  stands  for  the  following  time  series  process  model: 

''■'  2i;t 

x(t)  =  a^  +  a^t  +  -^  t(t-l)  +  (a^  +  a^t)Sin-^ 

+  (a^  +  agt)Cos  ^  +  a^Sin^it 

+  agCos^  . 

The  notation  described  in  this  section  will  be  used  throughout  the 
remainder  of  this  text,  as  well  as  in  Appendix  E.   In  addition,  all 
values  discussed  will  be  the  equivalent  single  smoothing  value,  regardless 
of  the  order  of  the  model  involved  (Section  4.4), 
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6.6  Wholesale  ChrysanChemura  Receipts  Time  Series 

All  tables  and  figures  are  contained  in  Appendix  E.   Hence,  in 
referring  to  a  specific  table  or  figure,  only  the  appropriate  table  or 
figure  number  will  be  given  to  avoid  unnecessary  repetition  of  the 
appendix  reference. 

The  wholesale  chrysanthemum  receipts  time  series  are  shown  in  Fig.  E-1. 
Note  the  growth  trend  of  the  curve.   This  tells  us  that  some  degree  poly- 
nomial should  probably  be  included  in  the  model  of  this  time  series.   The 
autocorrelation  analysis,  the  results  of  which  are  shown  in  Fig.  E-2, 
has  peaks  at  6  and  12  month  periods,  or  lags.   Thus  we  probably  need  to 
include  at  least  the  6-month  cyclic  function  pair,  since  the  highest 
autocorrelation  function  value  occurred  at  this  period.   Fig.  E-2  also 
shows  minor  peaks  at  21,  27,  and  34  months.   Due  to  the  limitation  in  the 
number  of  observations  available,  we  should  put  greater  weight  on  the 
lower  period  peaks,  if  they  appear  significant.   The  squared  nature  of 
the  peaks  indicates  that  harmonics  may  be  present. 

The  harmonic  periodogram  on  this  time  series  (Fig.  E-3)  exhibits 
peaks,  in  order  of  decreasing  magnitude,  at  6,  21,  8  and  10,  4,  26,  and 
17  months.   Again  we  reject  the  higher  periods  due  to  the  inadequacy  of 
the  data.   Notice  that  the  12-month  periodic  coefficient  is  of  a  very 
low  magnitude  relative  to  the  others.   Since  chrysanthemums  are  sold 
throughout  the  year,  we  would  expect  this  periodic  to  be  more  significant. 
The  6-month  periodic  has  a  sharp  peak,  confirming  the  previous  analysis. 
It  seems  that  the  12-month  peak  in  the  autocorrelation  function  plot  was 
merely  a  reflection  of  the  6-month  periodic  in  the  data.   The  harmonic 
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analysis 'becomes  unstable  as  it  nears  the  Nyquist  frequency  (Section  2.4) 
for  this  data,  as  we  would  expect.  At  a  period  of  35  months,  for  example, 
there  are  only  3  data  points  upon  which  to  base  the  value  of  R.   This  is 
obviously  not  an  adequate  representation.   Note,  however,  that  even  if  a 
periodic  does  exist  at  35  months,  the  period  is  long  enough  that  it  can 
probably  be  handled  adequately  by  the  adaptive  nature  of  the  general 
exponential  smoothing  process. 

The  various  models  tried  on  the  wholesale  chrysanthemum  time  series 
are  listed  in  Table  E-2.   In  analyzing  results  of  the  various  model- 
discount  factor  combinations,  let  us  consider  each  efficiency  criterion 
in  order.   The  model  efficiency  results  have  been  plotted  according  to 
three  general  classifications: 

(1)  Cyclic  models  and  combination  models  with  polynomials  of  degree 
one  or  zero; 

(2)  Strictly  polynomial  models; 

(3)  Combination  models  with  polynomials  of  degree  greater  than  one. 
This  method  of  classification  is  justified  by  the  similarity  of  the  plots 
for  models  within  the  same  model  group  above.  So,  we  shall  first  discuss 
each  efficiency  criterion  as  it  relates  to  each  model  classification.  We 
may  refer  to  the  model  classifications  above  as,  respectively; 

(1)  Cyclic  models, 

(2)  Polynomial  models, 

(3)  Combination  models. 

Table  E-2  shows  the  values  of  Ie],  mean  absolute  error,  and  a—   , 
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error  standard  deviation,  for  the  models  and  ^'s  tried.   All  of  the 
cyclic  models  exhibit  a  unique  minimum  \£\    for  a  .f  of  about  O.AO  or  0.45 
(equivalent  first  order^).   This  must  be  interpretted  only  as  an  approx- 
imate best  S   -value,  since  it  was  not  economically  feasible  to  cover  a 
large  number  of  values  of  the  discount  factor  due  to  computer  time 
requirements.   Nevertheless,  this  range  of  best-f's  is  considerably 
different  from  that  anticipated  by  Brown  (2,  179).   Brown  indicates  that 
a  £   less  than  0.7  will  rarely  be  necessary. 

The  order  of  efficiency  by  model  and  by  model  classification  is  also 
shown  in  Table  E-2.   The  best  of  the  cyclic  models  was  C12,  5,  PI.   It 
seems  significantly  better  than  the  other  cyclic  models  on  both  the  l£l 
and  Cg-  plots,  perhaps  due  to  its  relative  simplicity.   Though  adding 
amplitude  growth  to  the  12-month  periodic  in  this  model  did  not  improve 
it  over  the  original  model,  it  was  still  better  than  the  same  model  with 
amplitude  growth  on  the  6-month  periodic.   This  fact  seems  to  contradict 
the  results  of  the  harmonic  analysis  which  shows  the  12-month  periodic  as 
insignificant.   The  inclusion  of  other  periodic  terms  in  the  cyclic  models 
did  not  improve  the  |£l  and  <>"—  values.   One  model  was  run  with  only  a 
zero  degree  polynomial.   It  gave  high  values  for  the  efficiency  criteria, 
Indicating  that  the  inclusion  of  polynomial  terms  is  aiding  the  model 
efficiency.   Other  model  results  are  shown  in  Figures  E-4  and  E-5. 

Polynomial  models  of  the  third,  fourth,  and  fifth  degrees  were 
applied  to  the  wholesale  chrysanthemum  receipts  data.   Though  Figures  E-6 
and  E-7  show  very  similar  minimum  efficiency  criteria  for  all  three  models, 
we  may  note  that  the  fifth  degree  model  is  more  stable.   The  efficiency 
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values  for  this  higher  degree  polynomial  do  not  diverge  from  the  minimum 
value  as  rapidly  as  for  the  other  two  models.   This  is  a  desirable  charac- 
teristic since  the  consequences  of  choosing  a  bad  3  value  will  be  less 
grave  for  the  fifth  degree  model.  As  shown,  the  best  polynomial  model, 
P5,  produced  the  best  forecasts  for  fi  =  0.7.   It  is  doubtful  whether 
higher  degree  models  would  give  sufficiently  better  results  to  warrant 
the  expense  for  additional  computer  time.   It  is  not  surprising  that  a 
polynomial  of  such  a  high  degree  as  the  fifth  degree  is  necessary  to 
describe  the  time  series  when  we  consider  that  it  must  account  for  the 
cyclic  movements  within  the  data. 

Among  the  many  combination  models  tried  (Table  E-2)  C6,  8,  P5  gave 
the  best  efficiency  criteria  values.   Models  containing  higher  degree 
polynomials  were  generally  better,  as  shown  in  Figures  E-8  and  E-9.   The 
combination  of  trigonometric  fitting  functions  that  performed  best  was 
that  with  periods  of  6  and  8  months.   Recall  that  these  were  identified 
as  the  most  significant  periodics  in  the  harmonic  analysis.   The  assoc- 
iated 3  for  greatest  efficiency  with  this  model  was  0.6.   This  model, 
C6,  8,  P5  was  also  the  best  model  among  all  models  tried.   The  inclusion 
of  amplitude  growth  on  periodic  terms  was  detrimental  to  the  efficiency 
of  the  models  in  which  this  was  tried. 

6.7   Retail  Chrysanthemum  Receipts  Time  Series 

The  retail  monthly  chrysanthemum  receipts  time  series  plot  (Fig.  E-10) 
shows  some  growth  trend,  but  this  growth  is  not  nearly  as  explicit  as 
Chat  in  the  wholesale  receipts  time  series.   This  is  probably  due  to  the 
fact  that  market  expansion  is  not  as  extreme,  since  Manhattan  Floral 
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Retail  deals  mainly  in  Manhattan,  Kansas  alone.   Thus,  we  might  suspect 
that  a  polynomial  of  a  degree  less  than  the  fifth  will  perform  adequately 
for  this  data.   Again,  the  data  is  very  noisy,  since  we  are  using  dollar 
receipts  to  represent  sales  per  month.   We  pointed  out  earlier  that  these 
monthly  money  receipts  are  not  necessarily  for  sales  in  that  particular 
month. 

The  autocorrelation  analysis  on  the  retail  receipts  time  series 
(Fig.  E-11)  indicates  cyclic  behavior  for  periods  of  12,  24,  8,   and  6 
months.   These  peaks  on  the  autocorrelation  function  plot  are  not  clear. 
Their  squared  form  suggests  the  presence  of  further  harmonics  in  the  data. 
The  harmonic  analysis  on  this  time  series  (Fig.  E-12)  shows  peaks  at 
periods  of  6,  12,  4,  and  8  and  10  months  respectively,  in  order  of  decreas- 
ing magnitude.   Again,  this  plot  appears  to  degenerate  at  periods  greater 
than  35  months  due  to  the  inadequacy  of  the  data. 

The  models  tried  in  the  general  exponential  smoothing  of  this  data 
are  listed  in  Table  E-4.  The  best  performing  cyclic  model  was  C12,  6,  4, 
PI,  as  shown  in  Figures  E-13  and  E-14.   It  appears  from  these  cyclic  models 
that  all  three  periods  (12,  6,  4)  are  significant  enough  to  warrant  inclu- 
sion in  the  model.  Models  containing  periodic  terms  with  growing  amplitudes 
did  not  perform  as  well  as  the  same  models  without  amplitude  growth.  Most 
cyclic  models  performed  best  for  6  =  0.8.   The  best  cyclic  model,  however, 
performed  better  for  B  =  0.6. 

Again,  third,  fourth,  and  fifth  degree  polynomial  models  were  tried 
(Figures  E-15  and  E-16) .   Just  as  for  the  wholesale  receipts  time  series, 
the  fifth  degree  polynomial  model  gave  better  forecasts  than  the  other  two. 
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The  efficiency  curves,  hotrever,  are  very  close  for  the  fourth  and  fifth 
degree  polynomials.   For  these  models  the  best  B  was  0.8,  just  as  for  the 
cyclic  models. 

Only  one  combination  model  (012,  6,  P3)  was  tried  on  the  retail 
receipts  time  series.   This  model  performed  better  than  all  others  mini- 
mizing |e]  and  a—  for  B  =  0.8  (Figures  E-17  and  E-18) .   The  periodica 
contained  in  this  model  are  those  indicated  by  the  harmonic  analysis  as 
being  the  most  significant. 

6.8  Lily  Time  Series 

Figure  E-19  is  a  plot  of  the  potted  lily  sales  time  series  for  yearly 
demand  periods  of  14  days  each.   This  time  series  is  considerably  different 
from  the  previous  two.   Instead  of  an  obvious  trend  line,  we  see  a  devia- 
tion upward  from  the  zero  demand  line.   The  data  appears  to  be  strongly 
periodic  in  nature,  exhibiting  amplitude  growth  on  the  periodics  present. 

The  autocorrelation  analysis  of  the  lily  data  shows  peaks  at  14  and 
5  days.   The  peak  at  14  days  is  squared  and  unclear,  indicating  harmonics. 
The  harmonic  analysis  shows  periodics  for  7  and  14  days.   From  this  result 
we  may  suspect  that  the  5-day  periodic  in  the  autocorrelation  analysis 
actually  occurs  for  a  7-day  lag.   The  14-day  peak  in  the  harmonic  analysis 
is  also  squared.   The  7-day  peak,  on  the  other  hand,  is  very  sharp.   Thus, 
in  this  analysis,  as  well  as  in  the  autocorrelation  analysis,  the  14-day 
peak  may  actually  be  a  reflection  of  the  7-day  peak,  and  not  a  separate 
periodic. 

Table  E-6  contains  a  list  of  the  models  used  in  attempts  to  forecast 
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the  lily  time  series.   All  of  the  cyclic  models  (Figures  E-22  and  E-23) 
minimized  the  efficiency  criteria  at  ^  =  0.8,   This  best  ^  is  considerably 
higher  than  for  both  the  wholesale  and  retail  chrysanthemum  receipts  time 
series.   Perhaps  this  is  due  to  the  absence  of  excess  random  variation  in 
the  lily  data.   If  this  is  the  casual  factor,  the  lily  time  series  would 
not  have  to  respond  as  quickly  to  changes  in  the  actual  data;  that  is, 
the  lily  data  models  would  perform  as  well  with  higher  values  of  the  dis- 
count factor. 

All  of  the  cyclic  models  tried  performed  similarly.   Two  models,  C14, 
7,  PI,  G14,  7  and  C14,  7,  PI,  may  be  identified  as  being  slightly  more 
efficient  than  the  others.   Note  that  these  two  models  are  identical  except 
for  the  amplitude  growth  of  the  first  one.   Adding  amplitude  growth  does 
not  appear  to  affect  the  model  performance  one  way  or  the  other  in  the  long 
run.   This  is  somewhat  surprising  since  the  data  plot  (Fig.  E-19)  shows 
apparent  amplitude  growth  in  the  actual  data.   It  may  also  seem  surprising 
that,  in  nearly  all  cases,  a  first  degree  polynomial  performed  more 
efficiently  than  PO.   The  periodic  components  present  in  the  two  models 
cited  above  are  those  indicated  as  significant  by  both  the  autocorrelation 
and  harmonic  analyses. 

It  did  not  seem  feasible  to  run  any  models  composed  strictly  of  poly- 
nomial fitting  functions,  since  this  data  appears  to  be  mainly  cyclic  in 
nature.   Hence,  we  shall  move  on  to  a  discussion  of  combination  models. 

The  combination  model  results  are  tabulated  in  Table  E-6  and  plotted 
in  Figures  E-24  and  E-25.   Again,  the  two  periodic  terms,  14  and  7-day 
periods,  aid  the  model  performance.   Since  both  models  gave  similar  results, 
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the  addition  of  a  higher  degree  polynomial  had  very  little  affect  upon 
model  performance.   The  fact  that  we  need  polynomial  terms  at  all  is 
somewhat  surprising.   The  combination  model,  C14,  7,  P2,  was  slightly 
better  than  the  best  cyclic  model  noted  earlier.   It  produced  troughs  on 
the  efficiency  criteria  plots  at  B  =  0.9,  again  higher  than  the  wholesale 
and  retail  receipts  best  6  values. 

6.9  Poinsettia  Time  Series 

The  poinsettia  data  plot  (Fig.  E-26)  is  very  similar  to  the  lily 
data  plot.   Recall  that  the  demand  period  each  year  chosen  for  this  data 
is  35  days  in  length.   This  demand  period  always  terminates  on  the  same 
date  of  the  year,  December  25,  instead  of  the  same  day  of  the  week,  Easter 
Sunday  in  the  case  of  the  lily  time  series.   Because  of  this  fact,  week- 
ends will  occur  at  different  locations  within  the  demand  period  from  year 
to  year.   Since  deliveries  are  rarely  made  on  Sundays,  the  data  will 
probably  contain  a  7-day  periodic.   However,  this  periodic,  aw  well  as 
any  others,  will  be  obscured  due  to  the  "week-end  effect"  pointed  out 
above.   Hence,  we  may  expect  some  additional  noise  in  this  data,  not 
encountered  in  the  lily  time  series.   Just  as  for  the  lily  data,  we  should 
expect  the  best  poinsettia  time  series  model  to  contain  cyclics  with 
growing  amplitudes.   Recall,  however,  that  our  suspicions  were  not  entirely 
verified  by  the  experimental  results  of  the  lily  time  series  models. 

The  autocorrelation  function  plot  for  the  poinsettia  time  series 
(Fig.  E-27)  indicates  cyclics  at  periods  of  31,  33,  6,  and  28  days.   These 
periods  are  not  easily  explained  unless  we  consider  the  "week-end  effect" 
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upon  the  data.   With  this  In  mind,  we  may  probably  interpret  the  31  and 
33-day  peaks  as  reflections  of  the  expected  35-day  periodic.   Similarly, 
the  6-day  peak  should  probably  be  taken  as  an  indication  of  a  7-day,  or 
weekly,  cycle.   The  harmonic  analysis  (Fig.  E-28)  is  also  very  obscure. 
Nevertheless,  we  may  interpret  the  peaks  at  32  and  8  days  as  periodics  of 
35  and  7  days  respectively,  since  these  two  periodics  seem  to  be  the 
logical  choices  for  the  data. 

Table  E-8  lists  all  the  models  and  their  resulting  efficiencies  for 
the  poinsettia  data.   All  of  the  cyclic  models  used  gave  very  similar 
minimum  error  results  (Fig.  E-29  and  E-30) .   However,  some  models  were 
more  stable  than  others.   The  best  performing  cyclic  model  was  C35,  7, 
PO,  at  a  discount  factor  value  of  0.8.   Models  containing  amplitude  growth 
terms  were  generally  loss  stable  than  those  without  these  terms.   Still, 
they  seem  to  give  almost  equal  performances  at  6  =  0.8.   The  periodic 
functions  included  in  these  models  are  those  interpretted  as  significant 
from  the  autocorrelation  and  harmonic  analyses.   PO  gave  better  results 
than  PI  as  we  hypothesized  from  the  time  series  plot. 

Again,  no  polynomial  models  were  run.   The  most  efficient  combination 
model  used  was  C35,  7,  P2,  as  shown  in  Figures  3-31  and  E-32.   The  tendency 
for  this  time  series  was  toward  the  lower  degree  polynomials,  as  evidenced 
by  both  the  cyclic  and  combination  models.   C35,  7,  P2  gave  minimum  error 
at  6  =  0.9. 

The  best  model  overall  for  the  poinsettia  time  series  was  the  cyclic 
model  C35,  7,  PO.  We  must  keep  in  mind,  however,  that  all  models,  cyclic 
and  combination,  gave  very  similar  results. 
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Mean  pei-cent  error  was  computed  for  the  best  model  in  each  model 
classification  for  both  the  retail  and  wholesale  chrysanthemum  time  series. 
For  each  of  these  models  eZ  covers  a  range  of  6  values.   The  results  of 
these  computations  are  shov.-;i  in  Table  E-9  and  Figures  E-33  and  E-34,   This 
efficiency  crxtC'  .^^^rees  very  closely  with  the  results  shown  by  |c  j  and 

a—  .   For  both  time  series  the  combination  models  gave  the  lowest  £%  for 
3  values  of  0.6  (wholesale)  and  0.8 (retail).   The  mean  percent  error 
figures  are  relatively  high  due  to  the  very  high  noise  level  of  these 
data.   c%  could  not  be  computed  for  the  lily  and  poinsettia  data  since 
both  of  these  time  series  contain  zero  demand  observations. 

6.10  Notes  on  Harmonic  Analysis 

Throughout  the  data  analyses  above  it  appears  that  harmonic  analysis 
gives  a  clearer  identification  of  the  periods  in  the  data  than  does  auto- 
correlation analysis.   However,  for  extremely  noisy  data,  even  harmonic 
analysis  is  quite  obscure. 

The  analysis  is  performed  on  detrended  data.   This  detrended  data  is 
merely  the  residuals  from  a  first  degree  polynomial  regression  on  the 
original  time  series.   Since  the  object  of  detrending  is  to  remove  the 
trend  from  the  data,  some  time  series  may  require  polynomial  regression 
of  a  higher  degree  than  the  first.   Brown  (2,  72)  warns  that  detrending 
with  higher  degrees  is  likely  to  remove  some  of  the  cyclic  variation  in 
the  data.   This  is  undoubtedly  true  for  cyclics  with  very  long  periods, 
but  detrending  with  high  degree  polynomials  may  actually  aid  in  the 
identification  of  short  period  cyclics.   Such  was  the  case  in  this  research. 
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The  wholesale  receipts  data  was  detrended  with  second,  third, 
fourth,  and  fifth  degree  polynomial  regression.   The  residuals  from  each 
of  these  regressions  were  used  to  obtain  harmonic  analysis  periodograms 
(Figures  E-35,  E-35,  E-37,  E-38) .   The  higher  degree  detrending  seemed 
to  emphasize  the  major  cyclic  periods  and  de-emphasize  the  minor  ones. 
This  effect  is  seen  as  sharper  peaks  at  the  major  periods  and  a  leveling 
of  the  less  significant  period  peaks.   It  is  doubtful  that  a  polynomial 
regression  of  a  higher  degree  would  have  aided  the  investigation  further, 
since  the  fifth  degree  polynomial  fits  the  wholesale  data  trend  quite 
well. 


106 


7.   CONCLUSIONS 

From  the  discussions  in  Chapter  6  we  may  conclude  that  for  a  partic- 
ular time  series  and  a  particular  model  of  this  time  series,  there  is  a 
unique  value  of  the  discount  factor  that  maximizes  model  efficiency.   We 
may  define  maximum  efficiency  as  the  minimization  of  the  three  efficiency 
criteria  proposed  in  Chapter  6,  that  is,  standard  deviation  of  the  fore- 
cast, mean  absolute  error,  and  error  standard  deviation  about  the  mean 
error.   The  nature  of  the  data  and  form  of  the  model  are  major  factors 
contributing  to  the  value  of  minimum  error  discount  factor.   If  the  time 
series  is  very  noisy,  S  must  be  lower  to  allow  the  model  to  respond  more 
quickly,  and  vice-versa.   Moreover,  if  the  model  is  not  appropriate,  that 
is,  does  not  describe  the  time  series  process  well,  S  must  again  be  low 
to  allow  quick  model  response  to  changes  in  the  actual  data. 

For  the  two  chrysanthemum  receipts  time  series  studied,  the  most 
efficient  value  of  B  was  lower  than  expected.   For  the  wholesale  data 
cyclic  models  it  was  O.A.   This  fact  seems  to  contradict  the  findings 
of  Brown  (2,  116),  who  implied  that  6  should  never  be  less  than  0.7. 
These  two  chrysanthemum  time  series  contained  excessive  noise,  however. 
The  lily  and  poinsettia  time  series  models  performed  best  for  a  discount 
factor  value  of  0.8.   Since  these  time  series  contained  less  superimposed 
noise  than  those  for  chrysanthemums,  the  3  value  is  higher,  as  expected. 

The  inclusion  of  polynomial  terms  in  the  models  of  all  the  time 
series  studied  tended  to  cause  the  most  efficient  value  of  B  to  be  higher 
than  for  models  that  were  mainly  cyclic.   This  may  be  seen  in  the 
classifications  made  on  the  models  used,  that  is,  cyclic,  polynomial. 
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and  combination.   In  general,  the  minimum  error  B  values  were  smaller 
for  cyclic  models  than  for  polynomial  models,  and  smaller  for  polynomial 
models  than  for  combination  models.   It  may  be  that  combination  models 
were  more  appropriate  xi\   describing  the  time  series  encountered  in  this 
research.   Though  the  most  efficient  models  for  the  retail  receipts  and 
poinsettia  time  series  were  not  combination  models  but  cyclic  models, 
it  may  also  be  that  the  appropriate  combination  model  was  not  tried. 

The  model  classifications  used  in  this  research  were  justified  by 
the  error  criteria  plots.   Models  within  the  same  classification  produced 
curves  belonging  to  the  same  family  on  the  error  criteria  plots;  further- 
more, different  model  classifications  produced  different  families  of 
curves.   Models  within  the  same  classifications  also  produced  best  results 
for  a  unique  6  value.   These  phenomena  hold  for  a  particular  time  series. 

The  efficiency  criteria  plots  for  the  retail  and  wholesale  receipts 
time  series  were  quite  similar,  as  were  the  same  plots  for  the  lily  and 
poinsettia  time  series.   The  plots  of  these  two  pairs  of  time  series, 
however,  differed  from  each  other.   Recall  that  the  former  time  series 
pair  dealt  with  items  sold  throughout  the  year,  while  the  latter  pair 
dealt  with  seasonal  items.   The  inference  here  is  that  the  demand  charac- 
teristics of  the  product  studied  may  affect  the  proper  choice  of  the 
discount  factor.   By  demand  characteristics  we  mean  whether  the  product 
has  constant  demand,  seasonal  demand,  whether  it  is  a  luxury,  a  necessity, 
etc. 

Of  the  methods  used  to  determine  the  initial  coefficient  vector, 
one  method  was  slightly  superior.   This  method  used  the  coefficients 
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determined  by  the  detrending  analysis  as  initial  coefficients  for 
polynomial  fitting  functions,  and  zeroes  as  initial  coefficients  for 
all  other  fitting  functions.   The  use  of  other  initial  conditions, 
even  completely  erroneous  initial  conditions,  did  not  affect  the  best 
B  choice.   The  general  exponential  smoothing  process  adapted  well  to 
inaccurate  initial  condition  choices.   It  produced  efficiencies  for 
some  of  the  error  criteria  almost  equivalent  to  those  with  accurate 
initial  condition  choices  for  the  same  models.   The  main  effect  of  a 
bad  choice  of  the  initial  coefficient  vector  was  to  make  the  choice  of 
6  more  critical;  that  is,  the  trough  on  the  error  criteria  plots  was 
much  narrower  when  the  initial  conditions  were  improperly  chosen. 

The  periodic  analyses  used,  autocorrelation  analysis  and  harmonic 
analysis,  were  quite  reliable  in  detecting  periodics  in  the  data. 
Harmonic  analysis  gave  a  clearer  indication  of  what  these  periodics  were 
than  did  autocorrelation  analysis.   When  the  trend  movement  in  the  time 
series  was  apparently  of  a  higher  polynomial  degree  than  the  first, 
detrending  with  higher  degree  polynomial  regression  gave  sharper  and  more 
clearly  defined  peaks  at  the  major  periods  in  the  periodogram,  while 
leveling  off  the  less  significant  peaks. 

In  order  to  make  an  accurate  forecast  of  the  time  series  studied 
in  this  research,  it  is  felt  that  better  data  will  have  to  be  acquired. 
The  two  chrysanthemum  receipts  time  series  contain  excessive  noise,  while 
the  lily  and  poinsettia  time  series  do  not  cover  a  sufficient  number  of 
past  cycles.   It  should  also  be  noted  here  that  there  is  still  much  to 
be  learned  about  general  exponential  smoothing  before  this  technique  can 
be  most  efficiently  used  in  forecasting. 
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APPENDIX. A.   DERIVATION  OF  STATISTICAL  REGRESSION  EQUATION 

Let  us  suppose  that  we  have  a  series  of  n  discrete  observations 
y(t)  which  we  wish  to  describe  by  some  model.   We  assume  that  the  time 
series  y(t)  is  composed  of  some  process  plus  superimposed  random  varia- 
tion, (k-1) .   We  have  at  our  disposal  related  independent  variables 
X  (t)  which  we  may  use  to  estimate  y(t).   Thus  we  may  write 


y(t)  =  y(t)  +  £(t) 

=  6jj  +  Sj^x^Ct)  +  82X2^^)  +  .  .  . 

+  ViVi(^^  +  ^^'> 

where 

y(t)  =  the  value  of  the  original  time  series  at  time  t 

y(t)  =  our  estimate  of  y(t) 

x.(t)  =  the  i   independent  variable  at  time  t 

g  =  a  constant 
o 

3.  =  the  coefficient  of  x.(t) 
w^  i^ 

c(t)  =  the  error  of  our  estimate  at  time  t 
=  y(t)  -  y(t). 

We  may  write  this  equation  in  matrix  notation  as 


where 


Y  =  x6  +  e 


yd) 
y(2) 


y(n) 


=  an  n  X  1  column,  vector  of  time  series 
observations 
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1   x^(l) 


Vi^i) 


1   x^(2)  .  .  .  Xj^_^(2) 


1   x^(n)  .  .  .  Xj^_^(n) 


=  an  n  X  k  matrix  of  independent  variables 


"k-1 


a  k  X  1  column  vector  of  coefficients 


e(l) 
e(2) 


:(n) 


=  an  n  X  1  column  vector  of  fitting  errors 


We  wish  to  determine  the  vector  independent  variable  coefficients 
that  minimize  the  sum  of  squared  errors.   So,  let  us  first  determine 
this  squared  residual  sum,  S. 
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S  =  I   CD^  =  e'£ 

i=l 

where 

e'  =  the  transpose  of  e. 
Note  that  we  may  write  £  thusly: 

e  =  Y  -  XS 
So, 

S  =  (Y  -  X6)'  (Y  -  X3) 
=  (Y-  -e'X-)  (Y  -  XB) 

Multiplying  this  out,  we  find  that 

S  =  Y'Y  -  Y'XB  -  g'X'Y  +  B'X'X3 
=  Y'Y  -  2B'X'Y  +  B'X'XS 
=  Y'Y  -  2(Y'  X)B  +  B'X'XB  . 

If  we  find  the  derivative  of  S  with  respect  to  B,  and  set  the  result 
equal  to  Eero,  we  may  solve  for  B. 


II  =  -2X'Y  +  2X'X6  =  0 


X'Xb  =  X'Y 
Hence 

b  =  (X'X)"'''  X'Y 

is  our  least  squared  residuals  estimate  of 
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APPENDIX  B.   DERIVATION  OF  EXPONENTIAL  SMOOTHING  EQUATION 

Let  us  suppose  that  we  have  discrete  time  series  x(t)  which  is 
composed  of  some  tirae-model-dependent  process  plus  random  variation. 
This  may  be  written; 

x(t)  =  Kit)   +  e(t) 
where 

x(t)  =  the  time  series  observation  at  time  t 

^(t)  =  the  value  of  the  time-dependent  process  at  time  t 

e (t)  =  random  variation  at  time  t 

We  wish  to  approximate  f;(t)  by  a  linear  function  of  time,  x(t). 


x(t)  =  I  a.(t)  f.(t) 
i=l 
where 

x(t)  =  our  estimate  of  5(t) 

f .(t)  =  the  i   fitting  function  evaluated  for  time  t 

a,(t)  =  the  coefficient  of  the  i   fitting  function. 

The  fitting  functions  are  all  deterministic  functions  of  time. 
The  equation  above  may  be  written  in  matrix  notation  as: 


x{t)  =  a'(t)f(t) 


where 


a(t) 


a^(t) 
a^Ct) 


a„(t) 


an  n  X  1  column  vector  of  fitting 
function  coefficients 
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f(t)  = 


f  (t) 


an  n  X  1  column  vector  of  fitting 
functions 


Substituting  this  back  into  the  exquation  for  x(t),  and  again  using 
matrix  notation,  we  may  write: 

X  =  X  +  e 
=  a'^+  e 


where 


X  =  [x(l),  x(2),  ....  x(I)]  =  a  1  X  T  row  vector  of 

data  observations  up  to 
time  T. 

i  =  [x(l).  i(2),  .  .  .,  x(T)] 

=  a  1  X  T  row  vector  of  estimates  of  x(t)  given  by 
the  model 


^ 


f^(l)   f^(2)  .  .  .  f^{T) 


f2(l)    ^2^^'^ 


f^m 


f  (t)    f  (2)  .  .  .  f  (T) 

n       n  n 


an  n  X  T  matrix  of  fitting  functions 
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e  =  [e(l),  e(2),  .  .  .,  e(T)] 

=  a  1  X  T  row  vector  of  residuals 

=  X  -  X  . 

We  wish  to  determine  the  coefficient  vector,  a,  that  minimizes  the 

weighted  sum  of  squared  residuals.   In  exponential  smoothing  the  time 

It 
series  observation  k  periods  in  the  past  is  given  weight  =6  .   We  shall 

define  a  I  x  T  matrix  W  in  which  W  ,  is  the  square  root  of  the  weight 

given  the  observation  (T-i)  sampling  intervals  age,  and  all  off-diagonal 

elements  are  zero.   Now  we  may  express  the  weighted  sum  of  squared 

residuals  as 

S  =  eW(eW)' 
Since 

e  =  X  -  a'5^,  we  may  write 

S  =  eW(eW)'  =  (xW  -  a'^W)  (xW  -  a'S'W)'. 

The  vector  of  fitting  function  coefficients  that  minimizes  this  weighted 

sum  of  squared  residuals  may  be  determined  from  the  equation  (Brown,  2,  392), 


^=  0. 
3a 


So, 


and 


S  =  xWW'x'  -  a'S'WW'x' 

-xWW'S^-'a  +  a"^\n<''^'& 


3S 


0  =  -  xWW'^'  +  a'^W5^' 
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To  simplify  the  equation  above,  let  us  make  the  following  definition: 

=  an  n  X  n  symmetrical  matrix  of  weighted  fitting 
functions 
g  =  [-xm'^']'   =  an  n  X  1  data  vector. 

Now,  the  previous  equation  reduces  to 
g'  =  a'F 


a-  =  g'F"-"- 


where 


T  must  be  >  n  for  F   to  exist. 


This  vector  of  fitting  function  coefficients  minimizes  the  weighted  sum 
of  squared  residuals. 
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APPENDIX  C.   POLYNOMIAL  REGRESSION 

Suppose  we  have  a  discrete  ciine  series  of  T  observations  y(t)  to 
which  we  wish  to  fit  a  polynomial  of  degree  (k  -  1).   We  nay  express  a 
polynomial  of  degree  (k  -  1)  on  time  as: 

pCt)  =  1  +  t  +  t^  +  .  .  .  +  t''"-'- 

where 

t  =  time. 

This  polynomial,  evaluated  for  t  from  0  to  T,  may  be  compactly  expressed 
in  the  form  of  a  T  x  k  matrix  X, 
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,k-l 


k-1 


1   2   A  •  •  •  2 


1   T   T^  .  ■  •  t''"-'- 


where 


X. .  =  the  element  of  X  contained  in  row  i  and  column  i 


=  (i) 


J-1 


We  wish  to  describe  an  observation  y(t)  by 

k-i 

y(t)  =  E  S.t  +  c(t) 
i=0  ^ 

where 


y(t)  =  the  time  series  observation  made  at  time  t. 
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6.  =  the  coefficient  of  t  for  any  t. 
£:(t)  =  the  fitting  error  at  time  t. 

In  matrix  notation  this  equation  would  be 

Y  =  X3  +  E 
where 

Y  =  a  T  X  1  column  vector  of  time  series  observations. 
X  =  the  T  X  k  matrix  defined  earlier. 

6  =  a  k  X  1  column  vector  of  polynomial  term  coefficients. 
£  =  a  T  X  1  column  vector  of  fitting  errors. 

The  vector  S  must  be  chosen  such  that  the  sum  of  squared  residuals 
e(t)  is  minimized.   The  least  squares  estimate  of  B,  b,  is  given  by 
(Appendix  A) , 

b  =  (X'X)"''"  X'Y 

=  a  k  X  1  column  vector  of  least  squares  estimates  of  B. 

So,  having  determined  b,  we  may  write  the  resulting  regression  estimates 
as: 

Y  =  Xb 
where 

Y  =  a  T  X  1  column  vector  of  (k  -  1)  degree  polynomial 

regression  estimates  for  y(t). 
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APPENDIX  D.   THE  Z-TRA.NSFORM 

The  Z-transform  is  a  useful  technique  for  the  analysis  of  a  linear, 
discrete,  time-invariant  system  (Brown) .   By  discrete  we  mean  a  sampled 
data  system.   The  linearity  of  a  system  may  best  be  described  by  the 
following: 

For  input  =  x  (t) ,  then  output  =  y.(t). 
For  input  =  x,(t),  then  output  =  y„(t). 
Now,  when  input  =  a  x. (t)  +  b  x^Ct),  where  a  and  b  are  constants, 
the  system  is  linear  if 

output  =  a  y.(t)  +  b  y^Ct). 
A  system  is  said  to  be  time-invariant  if  input  x(t)  causes  output  y(t) 
regardless  of  the  value  of  t. 

If  we  let  f (n)  be  any  function  defined  for  discrete  intervals,  we 
may  define  its  Z-transform,  F(2),  as 

F(Z)  =  E   f(n)z"  . 
n=-~ 

This  transform  is  especially  useful  in  representing  convolutions. 

A  convolution  may  be  defined  as 

y(t)  =  I     x(n)  n(t-n) 

n=-to 

=  K(t)*h(t) 

where 

x(t)  and  h(t)  are  sequences  of  observations. 
In  the  case  of  system  response  analysis,  x(t)  may  be  the  input  signal. 
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and  h(t)  the  response  to  a  unit  Impulse  t  periods  after  its  input. 
Taking  the  Z-trans£orm  of  both  sides  of  the  convolution  equation,  we 
may  see  that  a  convolution  in  the  time  domain  is  replaced  by  a  multip- 
lication, a  much  simpler  mathematical  operation,  in  the  Z-domain. 

Y(Z)  =  X(Z)  •  H(Z) 

Y(Z)  =  the  Z-transform  of  y(t) 
X(Z)  =  the  Z-transform  of  x(t) 

H(Z)  =  the  Z-transform  of  the  impulse  response  h(t) 
=  the  transfer  function. 
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APPENDIX  F.     -      TABLES  AMD  FIGOHSS 

Each  timo   sorios  model   referred  to  In   this   appendix   5s   olasslflod   into 

one   of  tho    follov.'ino:   three   oatep;ories: 

(l)      Cyolioal  Models   -  modols  -.Thloh  contain  trlgouometrio   fitting   funo- 
tions   and   polynomials   of  def;ree   loss  than   tvro. 

(Z)      Polynotc.ial  "odols   -  models  which  contain  only  polynomial   fitting 
functions. 

(3)      Combination  Ii'odels   -  nodols  which  contain  trlf^onometric   fitting 
functions   and   polynomials   of  degree   greater   than  one. 

Model  components   are   indicated   according  to  the   following  notation: 


Cn 


P2.    Pk.    Pn.    °Pl.    P2.    Pic 


v^here  Pi,    n  =  integers. 

This   symbolic  notation   is    interpretted   as    follows: 


Cp, ,    pg,    Pj^  indicates  that   the  model  contains  both  a  sine   an 
fitting   function  for  each  of  the   integral   period 


and   cosine 
itting   function  for  each  of  the   integral   periods,    pi,    pj, 
and  Pv- 


Pjj  indicates  that  the  model   contains    fitting  functions   corresponding  to 
all  of    the   terms   of  an  n       dogree   polynomdal, 

''Pl>    ^2t    Pv;   indicates  the   inclusion   of  linear  amplitude   growth  on  each 
of   the    sine-cosine    pairs    in  the  m.odsl    "or  periods    of  pj,    ■po, 
and  pij.      These   terms  are   just  trigonometric  fitting  functions 
multiplied  by  time,    t. 

If  a  Eodel  does  not  contain  fitting  functions   of  the  C,    P,    or  G  type,    the   coi~ 

responding  symbol   (C,    P,    or  G)    is   omitted  in  the  notation  for  this  model. 


123 


TABLE  IDE;JTIFICATICN   CKAHT 
TABLE  HUldBER  DESCRIPTIQIT 


E-1  Monthly  Vniolesale   Chrysanthomun  Iteoeipts 

Actual   and    Da  trended    Time  Series. 

S-2  Wholesale   Chrysanthemum    Receipts    Tino   Series 

Models    Tried   and    fesulting   Efficiencies. 

E-3  Monthly   Retail   Chrysanthemum  Receipts  Actual 

and    Dotrended   Time   Series. 

E-4  Retail  Chrysanthemum  iteceipts   Time  Series 

Kodels   Tried  and   fiesultinj  Efficiencies. 

E-5  Lily  Daily  Sales  Actual  and  Detrended   Tine 

Series. 

S-6  Lily    Time  Series   I'odels    Tried   and    Pasultlng 

Efficiencies. 

E-7  Poinsettia  Daily  Sales  Actual  and  Detrended 

Time   Series. 

E-8  Poinsettia  Sales   Time  Series  yodels   Tried 

and  iiesulting  Efficiencies. 

E-9  Kean   Percent  Error  Values    for  Selected  i'odels 

of  "I'.'holesale  and   Eatail  Chrysanthemum  Tins 
Series . 
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FiaU:ffi  IDKHTIFICATICfl  CHART 

The  tabla  below  is  provided   as  a   reforenoe  to  all  of  the   figures  which 
follow. 

FIGURE  NI'MBER  DESC^IPTIOtl 

E-l  Vfholesalo  Wonthly  Chrysanthemum   feoeipts    Ti-e 

Series    Plot. 

E-2  Autooorrslation  Function   for  Tfholesale   Chrysan- 

thomum   ifeoeipts    Data. 

E-3  Periodogram  for  T/holesale   Chrysantheraun   Racaipts 

Data. 

E-4  I/ean  Absolute  Error  Versus   Discount  Factor  for 

T/holesale   Chrysanthemum  Data--Cyolio  Models. 

3-5  Standard  Deviation  of  Error  Versus   Discount 

Factor  for  V/holesale  Chrysanthemum  Data — 
Cyclic  l!odels. 

E-6  Koan  Absolute   Error  Versus    Discount  Factor  for 

Tfholesale   Chrysanthemum  Data — Polynomial  l^odals. 

E-7  Standard   Deviation  of  Error  Versus    Discount 

Factor  for  V/holesale   Chrysanthemum  Data — Poly- 
nomial Hodels. 

E-8  Kean  Absolute  Error  Versus   Discount  Factor   for 

Wholesale   Chrysanthemum  Data — Com.bination  Vodels. 

E-9  Standard   Deviation  of  Error  Versus    Discount 

Factor  for  V/holesalo   Chrysanthem.um  Data--Coa- 
bination  Models. 

E-10  Ifetail  Monthly  Chrysanthemum   tteceipts    Tina 

Series    Plot. 

E-11  Autocorrelation  Function   for   Retail   Chrysanthe- 

mum   fecoipts    Data. 

E-1?.  Periodogram  for   Retail   Chrysanthemum   Receipts 

Data. 

E-13  Mean  Absolute   Error  Versus   Discount  Factor  for 

Ifetail   Chrysanthemum  Data — Cyclic  i'odels. 

E-14  Standard   Deviation  of  Error  Versus    Discount 

Factor  for   totail   Chrysanthemum  Data — Cyclic 
Vodols. 
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FISURE  ia!!/BFJR  DaSCRIFTIO:! 


E-15  Mean  Absolute  Error  Versus   Discount  Factor  for 

Retail  Chrysanthemum  Data — PolyTionlal  Fodels. 

E-16  Standard  Deviation  of  Error  Versus   Discount 

Factor  for  Pjjtail  Chrysanthemum  Data — Poly- 
nomial l'*odels- 

E-17  Mean  Absolute   Error  Versus    Discount  Factor  for 

Retail   Chrysanthemum  Data — Combination  Vodel. 

E-18  ■  Standard   Deviation  of  Error  Versus   Discount 

Factor  for  Retail  Chrysanthemum  Data~Conibina- 
tion  I'odol. 

E-19  Lily  Daily  Sales   Tims  Series   P]ot. 

E-20  Autocorrelation  Function   for   Lily  Data. 

E-21  Periodogram   for  Lily  Data. 

E-22  Mean  Absolute   Error  Versus   Discount  Factor  for 

Ifetail   Lily  Data— Cyclic  r-Todels. 

E-23  Standard   Deviation  of  Error  Versus   Discount 

Factor  for   Retail   Lily  Data — Cyclic  l^odels. 

E-24  Mean  Absolute   Error  Versus    Discount  Factor  for 

Lily  Data — Combination  Models. 

E-25  Standard   Deviation  of  Error  Versus   Discount 

Factor  for  Lily  Data — Combination  J'odels. 

E-26  Poinsettias    Daily  Sales    Time   Series    Plot. 

S-27  Autocorrelation   Function   for   Poinsettia   Data. 

E-28  Periodogram  for  Poinsettia  Data. 

E-29  Mean  Absolute   Error  Versus   Discount  Factor  for 

Poinsettia  Data— Cyclic   I'odels. 

3-30  Standard   Deviation  of  Error  Versus   Discount 

Factor  for   Poinsettia  Data — Cyclic  Models. 

E-31  Mean  Absolute   Error  Versus    Discount  Factor  for 

Poinsettia  Data— Combination  Models. 

E-32  Standard   Deviation  of  Error  Versus   Discount 

Factor   for  Poinsettia  Data — Combination  Models. 
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FIOURE  NUMBER  D5SCRTPTI0N 

S-33  Wean  Percent  Error  Vsrsus   Discount  Factor  for 

Best  I'Jholssale   Chryaanthomum   Time   Series   J.'odela. 

E-34  Mean   Percent  Error  Versus    Discount  Factor  for 

Best   Retail   Chrysanthemum   Time   Series   i.'odals. 

E-35  Periodo,=;ram   for  TOolosale   Chrysanthemum  Data 

Usin,";  Second  Degree   Polynomial   .%?;rosslon 
Residuals. 

E-36  Periodogram  for  Wholesale   Chrysanthemum  Data 

Uslnp,   Third   Dajrea   Polynomial   Regression 
Residuals. 

E-37  PeriodO£;ram  for  TiTiolesale   Chrysanthemum  Data 

Usint;   Fourth  Degree    Polynomial    Regression 
Residuals. 

E-38  Periodosran  for  I'lTiolesale   Chrysanthemum.  Data 

Using  Fifth  Degree   Polynomial   Rsgrossion 
Residuals. 
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PERIOD  ACTUAL  DETRENDED 


I.                           27(4.C(;OCCCOO        .  23C.<;'t£21b't8 

27S.CC0CCC00  ,  i£e.ee'iuce23 

2S3.COCCCC0C  I23.t2t7e223 

lye.COOCCCCO  5.76904297 

373.C0Cr.CC0n  l'^7.  7  l  I2fi'i'.5 

261  .coccccon  iob.eG2S'i<;  i<; 

7,                           570.C0CCC00C  2  t  6  .  S«;'><3  A  72  7 

a.                   .1/.6.C000COCC  o9.s3aoe5<;4 

g.           566.0nCCCOOC  SSK.'.e0'46a7E 

IC.           A73.0CCCCC00  ii2.'.225<}2  le 

,1.           4C3.CCCnCCC0  -VC.63b0C')77 

12.           SAq.CCOCCCCO  ■    32.20712e91 


13.           £30.CC0CC0Cn  7C.215S117: 

^4.           7cO.C00CC0C0  1S7.15163513 

15.  ,  1  ICO.CCCCCCOO  ASA.  1337gQC£ 

16.  qji.COCCCCCO  2 3 J . C7 6 1  St  62 

17.  E4A.C0CCC0C0  -li;7.yeiA3CC5 
le.  1039.COCCCOOO  263.96069336 
IS.           521.0CCCCCCO  132.9C306091 

2c.                 ic.e.coccccco  206.eA5iQ9sa 

21.  1053.000(.COOO  142.75759766  ^ 

22.  ISA3.000CC00C  595.72973633 

23.  "  13eO.C0CCCC00  3e9.e7lU7SC0 
2i,.          1S09.CCCCCC00  475,61401367 

25.  IcSl.CCCCCOCO  574.556  64  063 

26.  1601.CCCCCCC0  4ei-49fc77930 

27.  2719.C00CCCC0  .  I  5 b t . 4  A  1  1 6 2 1  1  , 
15S4.C0CCC000  246.38305664  . 
1772.CCCCCCCO  S23. 32543945 

e76.CCCCC0C0  -415.73193359  ■ 

3i.          1 142.C0C0C0CO  -I92.7e977966 

32.  1404.C0CCCCC0  26.15234375 

33.  l64e.CnuCC000  227.C9A7112C 
3^.     ■         .     ie37.C0CCCC0C  173.03664998  . 

35.  1553.C0CCC000  SC.97924fcC5' 

36.  1177.0CC0C0OC  -373.07636914 
1361 .COOCCCOO  -232.13623047 
195C.CC0CC0C0  3l3.eOb66AC6 
1943.C00OCCOn  263. 74229102 
1337. COOCCCOO  -  J £S . 3 0932 6 1 7 

HI,                        1143. COOCCCOO  -622.36  7  16750 

H2.        1061 .CCCCCOCO  -727.42504883 

143b. COOCCCOO  -416.46266602 

IS'iA. COOCCCOO  -3AQ.S4U2a3?0 

i43d.aoccccon  -499.59790039 

1247.CC0CCtC0  -733.6So00bf 6 

;47.                         1  1  ?1  .CCOCCCOn  -652. 7136^305 

1291. COOCCCOO  -7?5. 77124023 

834.CC0CC000  -12/S.S2<;a5  7  42 

979. CCCCCOCO  -  1  1 7 3 . 6  Ufa  7 1 e 75 

51.  Iclfc.OCCCCCOO  -1 1 79.y44b«C0e 

52.  74H.COCCC00O  -  1  A v  1 . C 024 4  I  4 1 

TABLE  E-1 
Monthly  Wholesale    Chrysanthemum   Hsoelpts 
Actual  and   Detrended   Time  Series 


28 
29 
30 


37. 
38. 
39. 
40> 


43. 
'4  4. 
4E. 
46. 


48. 
49. 
50. 
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RIOD 

-ACTUAL 

DETHENDED 

53. 

1 '150.00000  COO 

-832.0598 1445 

HI. 

976.0000C00C 

-1349.  1  1 7675  78 

S5. 

ICOl  .CCCCCUOO 

-1367.  I75;j37  1  1 

.  56. 

1193.00  0  00000 

-1216.23339844 

57. 

10  2  4.0  0  000000 

-14  30.29077148 

58. 

Z028.COOCCOOO 

- 4 e9. 34 8 63281 

59. 

774 . OOCOCOOO 

-1 7Ce.40625CC0 

'  60. 

1  IIJI  .COOCCOOO 

-  1  422.40435547 

61. 

1356.0  0  0000UO 

-1270.52 197246 

62. 

2408. COOCCOOO 

-261.37958984 

63. 

1624  .COOCCOOO 

-1088.63720703 

64. 

1860. OOOCO COO 

-895.69531250 

6S. 

22 lO.OOOCCCOO 

-588. 75292969 

66. 

1701. COOCCOOO 

-1  140.81 079102 
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59.13159180 
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2404.00000000 

-523.92602539 
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2684. COOOC COO 

-2=6.98388672 

70. 
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553.9580078 1 

7J. 

2826. COOCCOOO 

-23 1 .09934998 

72. 

16S9. 00 000000 
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73. 

2C77.0000COOO 

-  1 066.21484375 

74. 

2658.00000000 

-528.272/0508 

75. 

4306. OOCOCOOO 

1 076.66943359 

76. 

3604. COOCCOOO 

33 1 .61 18164 1 

T'-. 

4 1 16. OOCOCOOO 

800.554 19922 

78. 

3102. COOCCOOO 

-256.503662 1 1 

79. 

2878.00000000 

-523.56127930 

80. 
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136  1  .38061523 
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3422  .COOCCOOO 

-65.67t74255 

82. 

4 196. COOCCOOO 
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84. 

3e63.0000CflOO 

46. 14988708 

85. 

3970. COOCCOOO 

310.09228516 

86. 

5S26. OOCOCOOO 

1823.03442383 

87. 

E575.0C0CCO00 

1828.97656250 

88. 

4956. QOOCCOOO 

1 166.918701 17 

89. 

3450. COOCCOOO 

-382. 13842773 

90. 

4620. OOCOCOOO 

744.60346680 

91. 

42S6 .OOCOCOOO 

337. 74584961 

92. 

4369.00000000 

407.687744 14 

93. 

5469. COOCCOOO 

1 4o4. 63012695 

94. 

6  117. COOOOOOO 

20b9. 57250977 

95. 

4546. OOCOCOOO 

455.Sl4a9i!5B 

96. 

34i:a. OOCOCOOO 

-645.54296875 

97. 

4  238.00000000 

6  1.40234375 

98. 

4200. COOCCOOO 

-1  9.C5625C00 

99. 

4C3 1 . OCOCCOOO 

-23 1 . 7 1484375 

100. 

534  I . CCCCCCOO 

1 0  35.22656250 

10  1. 

3712.00  0  00000 

-636 .828 125C0 

102. 

369  7.00OCCOC0 

-694. 8867 1 875 

103. 

4673. OCOCCOOO 

238.0546  8  750 

104. 

3557.00000000 

-92  1  .00390625 

105. 
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PERIOD  ACTUAL 
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fic-.  r;-3s 

Using   Third   Ds.-jroo    i-'olyncr.-.l.tl    p3!;ressio?i   R'isidianls 
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Of  the   proj^rirr.s    listed   in    i'vis    Jippandi:-:,    Al  i'CCET,    PCVSTSC,    and   EXPOS"0 
ai-o    rcsvisions   o"  progr?j?3   originally  v/ritton   by    fhj-nioni   Miller    (7,    65).      Tlie 
Con;putsr   prof;;rr.ns   usod  vere: 

(1)      M-ICDSri 

(£)      l'0',V3Pi;C 

(3)  AUTCPLT 

(4)  PLCTOOl 

(5)  expos:  :o 

(fl)      1-ATREGH 
AUT0D3T  p2rX'orm5   E  first  do;-;re3   polynoirdal   rsgressicn   on    lh3   orii^iinal   time 
serios.      Tiio    resultiuj   residuals    fonr.  t:is  dstroAded   tiina   series    f rem  v.'hioh 
tho   pro»;rarr!  ootr.putos   an   autonorrolation  function   for  spsoified    lag  valres. 
PCTSPSO  uses   tho  dotrended   timi   sories    from.  Al'TCDBT  to   perform  a  har.'Tionia 
analysis    for  spooifitid   poriods.      It   computss    ths   coefficients    of  har?-ionios 
in  the   arithrietio   trii^onoir.etrio   sequence   explained   in  Section    £.4,      The   re- 
sulting ooeffiCiient  values  nay  be  used  to   construct   a   poriodogran  for    the 
tice   ssriss.      AVTOPLT  plots   the   orlf^inal   and  dstronded   tire    series    froa 
AUTODET.      PLCTOOl   produces    autocori'elatjo:!  functi.on  and   periodor.ram  values^ 
ooDiputod    in  Al'TODST  and   PO..0PEC,    respectively. 

EXFCSrO  perfor.-iS  the   jjeneral  exponential  STiOothing  process   on  the   orig- 
inal ti;ii.?   series.      At  each  samplinf;;    jxitprral  a    forecast   is   T^ade   of  the  next 
tire   series  value.      Each  fcrncast   is    based  on   all  data  vo  to  that  particular 
point   in  time.      In   addition   tc  the    tin-.s   series    itself,    tha    following  data  must 
be   fed   into   this    proc;ran: 

(1)      !hinb^?  r  of  observations   and  number   of  .T^odel  fitting  functions. 
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(?)      Initial  ooaffioiant  vaotor. 

(3)  Input   fotT.iat  for  trfjisition  matrix. 

(4)  Transition  ratrix. 

(0)      Irput  fonnat    for   time   series. 

(g)      Valu'js    of  discount    factor  to   be   tried. 

For   a  mora   ooipprshoisivo   explanation   of  AIJTODET,    POV,'S-EC,    nnd  EXPCSVO 
refer   to    tha  master's  thosis   by  Millar,    (7,    C5) .      liATRE'GR  is    j'jst  a   general 
regression  pro^rsjn*      T?i9  niatrix  of  independent  variables  nay  be   of  any  desir- 
ab}e    form.      This    proj^ran  outputs   &s\   analysis   of  variance  table   and   a   plot 
containint^   the   actual,    ostirriated,    and   residual  time  series. 
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cor.'PVTr'R  noo;i!.r  ids'itifioatio;'  cvi-.ri 

Tho    list  ^3^ovf  is   provid-3d   cs    r.    rsforc^ioe  tc  all   of   tha   oon'.p-jter   prosrairs 
which   folloTx: 

F_l  AUTODiST  -  detronds    the   orlsinsl   timg    sgriss 

and   ooiTiputes   the   autocorrelation 
function, 

F-2  POVfSPSC  -   ooKiputos   psriodio  term  ooaffioionts 

for  an  arithmotio  sequence   of  periods. 

F_3  AUTOPLT  -    Plots   original   and   detrsnded   tiyre 

sories    from  ATJTODST. 

F-4  PLOTOL'l   -   plots   o'.itooorralation  function   from 

ArTODET  and   pariodogram  from  PC'SPSC. 

F-5  EXPCSrO  -  Exponentially  smoothos    the  original 

time   series    for  any  model  and  dis- 
count   factor. 

F-6  MATffiGR  -   pprfonns    a  inultlple    regression   for   in- 

dependent variables    of  any  desired 
form. 
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11    Tfj-r../,  r  ( ;;  TH    thI'j.    ts 
1    hf,-:r.v.  r  ( i.T. :■.).-,) 
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_  Cir-'END  1CNX(  iGC  )  .  1  (  I  CC  )  ,C  (  1  CC  )  .-J  (  ICC  )  _  _ 

I   fr.rim{]2)  •' 

i'     I-  L.  r !  ■-■  (.  T  (  :■  F  2  C  .  ?  .  2  I  1  C  ) 
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tit     CGNrlNUli 
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PSDGiV.I.'  NO.   F- 
POr.'SPSC 
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DIMENSION     XCiCC  )  .  CUT  (  101  )  .YFH(  I  n  .  ANG(g)  _. 

1  FG"!VAT  (212) 

Z  FCr-^AT  (ZF2C.e) 

3  FC4.VAT  (  £5>  .  37hACTL."- (  1  )  .  C  ■;  T  k:^.NC£0  (  2  )  CATA  VS.   TIVE/) 

4  FCHN'AT  (F5.C.2F2C  .e  ) 

5  FC«r-'AT  (SI-     T  I^E  .  lOX  .eFACTCAL  .  1'4X  ,  3t-CElBEiNiDtD/) 

WFAD     (l.l)        S,NL  ^        

HEAD     (1,2)      (  X  (  I+2*N  )  ,X(  i  *N  )  ,  !-  1  ,N)         ' " 

CO     ICC         1= I .N 
:OX(I)     =     I  ___  _  _  1_ 

ViRITE     (3.5)  "  ~ 

URITE     (3.4)      ( X(  I )  .X (  I  ^^)  .x(  I  t2*K)  .  1  =  I  .N) 
VlR  t  TE     (3,3) 

CALL     PLOT     (  I  ,X,N,3.NL ,0 )  ~ 

ST(JP 

END 


PRCarii  IJ6..  F-3J 
AUTOPLT 
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c   Fi  )   PL'JT    cr    itTC    rorFTi.  •    ■■ .' 

c   h-   ■=   L'vr-:,      .■.   =    1    i    .M.  r'j  cci:f. 

C     i\    -      PLLT      TK       JC.■.:=^^';C      '"ATA 

C   LAG    =    vatr.     \=    l,-;    =    2'3    =    -i        - —         -•    -•   - 

IM  vl:J  ti:\    ,'(  1  -G  )  .     J  (  ^■.  C  )  ,L,i,  I(  1  C  1  )  .  v.-  •  (  1  1  )  .  AKG  (-7) 
1      f  0  .'•    AT      (  I  J  ) 
■"       2    FC'7;'..T     (r  1  c  .C,F  2C  .2  )  --..--.     

.J    fun,' AT    (  31= -c  .  y  .r  I  u.c  ) 
'(    r  (J '■.■.■  A  T    ( r  1 0  . 0  ■  3  f .-  0  •  7 ) 

5  rOii'  AT     (  J).  ,  !.IL  1  G?  .  1  C/ .  :  Oh  At  TC    CC^^FE/)  ' 

6  Fr,).'»r    !  I  (.t    F'<  fcLiiivCY  ,  1  (.' < ,  ;fa    ,  it'.<  .^j-.^j    .it.x.iHW   /) 
';«:^n    (  i  ,  i )    r^ 

"    '      -  >='C\n    (1,3)     (.M  [).!'(  I  ■••-),  1  =  1  .N  )         -         

>■  ■<  1  T  =     <  .3  ,  -3  ) 

i-:  ■;  1  T  S     (3,2)      (  P  (  1  )  ,  |J  (  I  t  N  )  ,  I  =  1  ,  i\  ) 

Call    '""Lcr    (  l  .P  ,:s,I',  c  ,  c  )    "■  

i'C4C     (  1  ,  I  )      \ 

'(i4n      (1,J)      (  A  (  1  fN  )  ,  A(  I  +  2i.'0  ,  -.(  ltjt\  )  ,  A(  I  )  ,  1  =  1  ,M) 
'  .■.'Jl  T=     {  3,  J  )                                                                                                                               ' 
VM-iiri      (3,>»)       (    •  (    I  )  .A  (  I  +\  )  ,  '<  (   1  +:?<  N  )  ,  A(   H  3^..N  )   ,  !  =  I  ,N  ) 
CAUL     PLUT     (  2,  A,,s,  4  , /-J  ,0) 
STOP  — ■■- 


PBOGFat.'  HO,    F-4 
PLOTOOl ' 
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[F(i-:n  )        1  3.  1 .1,  1  4 
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i)    (cun.Ai 

J)        FMT 
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:hk( r ) , 


1=1 ,\) 


J) . 


1=1 .N) 


J) . 


;:-^i  r!-;(  3,  5)  '    ^ 

DO     15        I=1,N 
LJ_''i^  [  Tr:  (  3  ,  6  )         (  T  ■"  (  I  , 

r^FTViJTiVT)       xFvT 

C    .iF.AC     CJ     THr.    Q gSEKvA  T  ICNS 
;JEAD(  I  ,  xf;:t  )        (X(l), 

rcu 


1  =  1  .,JD) 


PI 


ii: 


=     3.141: 
■.n(  1  ,'4  ) 


ah;-,  fa 


IF(  4f;CT4)         1?  ,  12,ri 
C        r.AV;-.    TM1-:     INITIAL     VALUCS    CF     TH.f    CUiNSTAryTS 
^ ^^K ^lli  _L'i' ^  '  * '-     V*l-iJ,-S     CF      rnu    PITTING    FUNCTIONS 

II   no   16     1  =  1  rs^ 

A(  I  ,  1  )      =     A0(  I  ,  1  ) 
^l  c     Cl(I)=Cu[(I) 

~calC 


. < ,C1 , 


■*AY    (.-ao.MOij", 
CALLMATINVCF ,N) 

CALL^Vf-:C(F.  \C  ,  N  ,  H  ,  jet  T  \  ) 

C  4 1  <  «  vF  o  Y"    THC    V    :'Af'Tfy 
CALLV  -A^!!  ■<(  FK  .F  ,rj) 
CALL     FOMC-.T     ( -JiND  ,CI  ,  A  1  ,  > 


.CI' 


.  ry  ,ci3cTA,f"k  , 


AG,  Al  ,  Ae>£TA,CO' 


i  i'idJX  ,^■■),  T'- 


UC    SO        II     =     l..\0 

rJLI   (II)        =        II 

PLT  (  IJ^^^lO  )     =     x(ll) 

"pltTi  I  flTS^jT)    =~-friTTri~) 
"LT(  1  It:^i.^;n)    =   xdii    -   fcsti 

50    I'l.T  (  I  I     t     .i*.NU  )     =     0.0 


.F .H.FCST) 


I) 


CALL     tJLOT 
KNN  =  \NJNt  1 
CO  J C     10 

2  CO'\tTJlJF 
ZC,     TO     t) 

3  'THi' 


(N-. 


.  PLT.MC  ,5,--IL  ,0  ) 


LSil 


PROGtoL'  SO.    F-5 
EXPOSI'O 
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~.U.:!''<C'UTI  ^•^     7AY     (^0  ,\TCL  .  -JtA,",  .Tl  .   'N  .  CH<  ,  r  V  i  CC-T  A  ,  Fi<  ,  F  ,  AC  .  A  1  .  ALi-^ 


ICuI  ) 


bi  ■•K-Nnic-;   r:<(  10  , 1  u)  ,  Ti  \  (  1  c  ,  iG) 

D I  ,-.::.MS  r  DJ    \  (  10 ,  I )  ,  1  (  1  ,  1  ? )  .  :  u  n' .  1  n  )  ,  F  {  I  0 .  I II )  .  C.J  I  ( 1  0  ) 

L;I''.i:;>:3ICni    <(  iici )  .ci  (  l  o  )  ,  ^  !  <  i  j  .  i  )  .h  (  i  o  ,  i  )  ,  ^c  t  i  g  ,  i  )  ,  r."  (  i  o  .  i  o  ) 


rjiMsrjsio  1  Ai=(  Id.  1 ) 

COU  JI.!;    P<=CIStlJ\     F.FK.CI.Cll 
I     FO^.v.AT     (•     F-v.ATt!l  X  •  ) 

AP  =  ( 1  ./rN) 


c   'i^^^rI^:■i    h:..j^    r:-t    F    -.-iTwix    is   '-fsc ALCULATtc    .virn    a   ne^,    sSta 

67    06-:  rA  =  (  A^^-:  TA  )  i«AP 
C^  ■!':  IA.-I  r  I  ALi  .i'l     THE     VECfG.?     CF     FITTINo     hUNCTIO^S     ANO     CGMSTAMTS 

C[(I)=CQI(I) 
42     AA( [ , 1 )=AC( I . I ) 


C    CUKPUT.:      THci    F     MATRIX    h\Jf!     ONE     MGOt:ir~AND~ffNT~STfA~ 
r)i;;i  I  =  1  ,W 
DL|- J=  1  ,i-i 
f^{  I  .  -J)=C.  0 

li   r  (  I ,  J  )  =  0 .  c 

RCFA= 1 .0 


CNr=  1 . 

;-3    DO  !0  [  =  1  .N 


30     J  (   1  .  t   ) =A (   I  .  1  )  7- 

eo'to  1-1  ,N 
:•  =  I 

ci<  I  ■  J)  =  M  I  .  i)fd(  1 .  j)*(-^i:ta) 

10    CQMTI  ^ll;E 


DD43I=l .\ 
DJA-lj  =  r/ ,  \ 


43    CI  1  <  I  .J  )=C  1  (  1  ,  J  )*l.STA 
I  F  (  K  -  I  )  o  J  .  1  J  ,  o  6 
_£*^*j;;fCHjx:<  JI^SC     CC -^V^IPOEJ^CE     OF     THC     F     t.'«T<l:<     AT     liACH     SOTH     ITC'JATIUN 
C«*i.*i.ThIS      «:uri.C     CmTzckS     CnTi'l     Pr^l^i  v3s~ArT"~cUET'HNT^^vhlcT~"iS^0T~T"f  m^ 
C*t.l-.*<-I  ZFi) 

c(i    iF(  c  •:  r-50  .  1  1 1 .4  4  . 1  a 


4  4     C  N  T      =     0  .  0 

DO    30        1=1 . M 

f  =  I 


OOiC  J=:',  \' 
u-:n:'  =  i-(  1  ,  J) 
A.MiJ.v^Cl  (I.J) 


•^A  r  In^ANU.'/'JdMI-r 

1K(  R  •.  T  lu- .  CCCCO  1  )  OO.KC  .  16 

-ctt**iTMr.   F    '/Ar-;i/    vlZijt    Tht    cc<vFq,:,7-.c.-:    c^^iteria    .vhfn    rms  cccu>!s    ri-e 

C  *  *  '■ '  »    F '■  .^  T  ^  1  N'     IS     CI,'- PL- TEC     .3y     CYl,•^;:rvY     M-,0     Ti-'E     P^CG'JA:.,     ChEC<S 

C  »  4  *  *  t  T  O     3  h  E     I  F"~fT,T:~S  A  NGE  '  uF~";3lrf  A^'Tl  .<  ;.~rE  tKr^a~irfSlrD 

0U.3  7I=^.^J 

i_=I-]  EXrcS).-n   (cont'd.) 
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D0a7J=l ,u 
3  7   r  I  I  .  J  )  =  F  (  J  ■  I  ) 

:■  •-•  I  T ;-:  I  3  ,  1  ) 

2ti     i-   <lT.t  (  ■!.  i  )         (f-(I,J),      J=l,\) 
GO     TP     'f'J 


Cttt**CKl  Ti=^  f*     IS     -Kir     I-LT     ^NI)     tJQl^KAL     lTii^AriG\     CC^TI:^lUeS 

10   ~(^50I  =  ^,;^! 


.v=  f  " 

1>050  J  =  --!.\ 

FK(  t  .  J  )=FK(  I  ,  J  )  +C  I  1  (  I  ,  .1) 


SO     f-  (  I  ,  J  )=h  (  I  .  J)  to  1  (  I  .  J) 
Cf  +  ti*^C./     unCMi;     Tr-.f:     vrcrG<     r.K     FITTING    FUNCTIONS     USING     THE    L     MATRIX 
0031 1=1 ,N 


.3  1     A  f  (  t  ,  1)  =  0  .  0 
00  >0I  =  1  .  ^ 
OIJ50  J=l  ,■•! 

001001=  1  .N 

AA(  1  ,  1  )  =  ^F(  I  .  1  ) 


ion     A( 1 , 1 )=AF( 1,1) 
C?:»**»THIS     IS     \     •JOtrli-vL-     FOv     CD  JVi;  .iTlNG     THE     VICTOR     Cr     FITTING    FUNCTT'. 
CJ^S  *«*F;j.J     T     TG     THF     v^rCTGJ     -;  V  VLU'^  T 10    KCH     -T 
lJoTJTi^T.T-. 

CK  =  Chi;(  I  ) 

rF(  CK-  l.)63,i.i.t3 


S'j     ^(  I  .  n=-i(  I  .  1  ) 

63   Ci' s:ti<^iui; 
C    SAVC    F(l)     FL)-1    fg:-!i:cast 


IF  (  K-  1  ju^.ai',  >j;d 

2i'    ()G6;!I  =1  .N 

^l(I.l)=^^(I.l) 


"CCTr  I  ltuT2 
F-rTA  =  [)L',£rA»Bi:TA 

c-jt=cmti-i  . 


i^-K  (■  1 
GO     Ttl    99 
Ij    CU'ITI-JUE 


EXFOSHO    (cont'd.) 
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r.Ui.HQUT  l-iZ:-  \  T  riv{C.'<) 
"T I  'iu^rn Tc;  J  "cTi  i; ,  1  vT,'^  ( Vi^)  ~,V  (  i .. 
O'lUiiL-    M'.-c  I  s  1  jt,   r.F-'<  ,ri  ,C!  1 , 

ALL=C( 1,1) 

OCl     /w. .',  c     J=l  ,\.vi 


•1  '■  4 (,    V  (  J  )  =  cJlTTTR T7 ALL" 

V(N)=I ./ALL 

DO     ')Ci<10     I  -1  ,N(.'l 

IT'TS  iTT 

CI  >l  1  =C(  I:j  I  ,  1  ) 
iJU    y:i-jO    j  =  1  ,:.r-i 


a  j -J 0  c  ( I ,  J )  =  c !  r^i"r.rnrT~:^^~c  i"PTi'* vi  xr 

■JT5n    C(  I  ,M)  =  (-1  .C)4CIPl  1  tV(N) 
DO     f)i,6  0     J  =  l  ,N 


>co  c(^rTri"vrTr 


_EXPOSiro   Co_ontld..) 
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su:-.>^i;ijT  I  i.7    V?- Aril  y.(FK  ,r  ,.\) 
7ri'7^fjTriT,N~F (iimoi  .i^Ki id; TO) 

c;l  vF^i^  !0  JVJ  {  1  c.  ic  ) 
oi'E-i'".  1  UNVl  (  I  >: ,  1  u) 

m  re -ca  r  ( ///I  o.< .  JvM   THr    7A.M>Kc:;   or   thi   cct£Fr  ic  i:z\rs    i; 

oil  i^yiTT^'i 

DUI4 J=l .N 

I  -'»    V  s  t  I ,  J  )  =  0  .  0 


oci  1  (-1  ,^i 


vj  ( ;; ,  I  )=■--(  K  ,  J  )  *FK(  J,  I  )  +  V2  (K  ,  I  ) 
OtJl  3  1  =  1  ,N 
UIH  3J  =  1  ,\ 
Vl  "(TTTr57r."l) 
'JO  121  =  I  .N 

Df)iaK=i  ,N 


0U12J=1,M  ' 

V I  ( K  ,  I  )  =  1/ ;:  (  < .  J  )  *F  (  J .  I  )  +  V  H  <  .  I  ) 

'.•;.<  I  TF  (  3,  7  /■?  ) 


OCI12  2v<=l  ,  ■•) 

if.t^lTlK  3.  7d  J)  V  1  (  K.K) 

RETUWJ 


liMO 


„  JXPOSJAQ_i.cfint_'  d , 
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H  V  i;c  ( r ',  A  c .  N .  H ,  11  .-■  .1  r  4  ) 
T-Ti 0 . n:T taq rr(r.T jtm' ( t o7 n  vi^f a u v r s 

;CIUI0  w    F.FK.Cl  ,ci  1 


"■.u  j-'iiij  r  I  :i 
I    Fc^-^riA  r    (  • 

7^7     HC--    'AT  (  /,'/  1  X,  lf;h      T. 

4'i't   rtr^MAT    cm    ,  i  ooi  3. 5) 

OSq    FC-l".Ar( ///iOX.6H     G£TA=.!^9.f>) 
l)03il  =t  ,'-1 

(1" 


OIVcNSIC-J  SvcoTHInG  1' 


vcCTi: 


0 


:  I   .N 

a  .-I 


'j13 

O  i?  'J 


H(  1 .  n  = 

0Lyj43J  = 

^R't"rr)"^ii~rrnT~r"rrrJT*ATr(~j"m~ 

CCNTINUE 
*V.^ITi:     TMi    CFFiECTIVc     oETA 


«T-MfF(3,l)        N 

vr^l  Till  1  ,;iS9)Gt;i:r4 

-.'■lITt:    F     l;j./E:^SF 


3.3  J 
~C"~ri7T 


DIJ  3.i  3t  =  1  .  M 
v.  M  T!i(    !,  i':4  )  (i 
JT 


(  I  ,  J  )  .  J  =  I  .  N  ) 


■) 


-;  I  T';;(  3.  7ci>(  n(  I  , 

ETUtJW 
NO 


EXP03V0   (cont'd.) 
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SU'Ji^dUTlN:!    FC>!CSr     (  N  ,  NO  i  C  I  .  »  1  .  ■<  .  *nfi  T  A  ,  wi-l  ,  Ti-!  .  F  .  H  ,  FC5T  ) 

~o rMfZ N'si "d'T *F c s"r c'To 5 o"y , >:TT{ so bvirc  i av i o y "rc c f rroT .  vr(  ro T i )' tat i o ,t r 

DI>;EN=;I0  JX(  ZOO)  ,H(  10.  1  )  .Cl  (  10  1  .  T:-'(  1  0  .  10  )  .Ci  F(  10  )  ,  CIL(  10  ) 

ui'.'EN'iin  ;   T  T-- ( ic,  10) .    PCTt:«(20o) 

~oou>rDErT=  i"^:i2ricK^rVF<  ."cr.cTiT  '         ' 

1  FO-;:iAr  (•  sto  u^v  FG^^^CAGrs  i -s  '.hi  a. it/) 

2  FO-e.'-'\T   (•  rSAU  E-Ji^O'*  IS  '.rli'.'i// 


3 

Ki'^.VAT     (•     MEAN     AL'SQUOTil    ■i-lKCK     IS     '.Fia.o/) 

4 

FUl'»MAT     (•      VA<14NCe     OF     THE     E.fnC.J     1-5     •.F15.6/) 

b 

FOi^f'.AT     (■      STD    DK«     '.t'CLT     NCAN     CT     ThJ     EfiiJC-i     IS     '.FlS.e/) 

6 

7 

FWHAT     (•      S  ro    UCV     AfiUUT     ii-J^a     (If     TH-     criUd'^     IS     *.F12.0/) 
FORMAT     <•     .-F.  \N     PFliCFNT     E?.-'C:<     '.  ^,      ".ri2.S/) 
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FU  ?;.  AT( ///I  <  .  3U-     THI-;     VAiMAf.C^"     GF     THE    FO:^EC  ASTS=  ,  F  I  3  .  5/ 1  h  1  ) 

107 

|-ov,iA  r(  l-;.F  jc.  7  .rao.a  ) 

r 

lOU 

FQt^fVArt  IX,  ntl     NEU     FQ'-JECAST            ) 

101 

FQ ?MA r( IX , 1 IH    rCDEL     NO . - . I ? , 20X . 6H     aETA=,F4.3) 

lya 

FO-»r.A  r(  1 /,i:Oh<">*TH£    Efii-iU)     IS***                         ) 

urn 

FO:^MAT(ajh        PCMCD                              Oli  S^:'<  V  A  T  I  CN                   FCHECAST 
1               CUM.     Fiii^GR              nCT.     Ehfi-iCl!      ) 

EUhJOfJ 

c 

lO'i 
ICn 

FC'"<,'"  A  T  (  1 .1  ,  F  I  J  .  I  ,  F  I'j  .  2  .  F  1  ■i  .  2  ,  P  1  -J  .  2  .  F  1  o  .  2  ) 
Ftl-!liAT(24H     SUV     CF     -SOUA'JE     OF     d  R,<a-<=  .  F  1  5  .  J  ) 

pmNr    CUT    THE    MnaEL    nc.    aud   heta    fck    this    F    -UTWIX 

WrJtTt:(  .1.101  ).V|J,  AMET  A 

c 

ViAKi 

-:     THE     Tf(ANSPGSF     CF     THE     Ti^AlS  I  T  I  Hi  1     MATRIX 
SUM=O.C 

su;';a  =  o.  o 

5Ui'!X  =  0.0 
SU^'lMtl    =     0.0 

SUVPCT     =0.0 

5EDSC     =0.0                                                                                                    .                                                    . 

Oa42'lI=l,\ 

I 

42M 

0G42'>J=  1  .N 

TTM(  I  ,  J)=Tf;{  J,l  ) 

c 

r-  AM 

-:     THE     FHTiEC^ST 

Wrt I rE( T.O  ?  1  ) 
'Otl  l'HK=  1  .  -JO 
t  =  0.0 

LIU'l  !  t-  1  ,M 

ZM  =  C [ (  I  )iA 1  (  I  ,  1  ) 

Tl\ 

z=-r  +  ;;.j 

Fcar(  K  )r 
c  calculat:-  th=   z^.'A}'i 

OHK)     =  .<(a)  -  I 

;5cTF'vTi.r)=^'^CeRCK:i)/Xi^^>  ;^  1  JO  .~r 

C     CALCJLAT.:      TH-      Cir-.L.l.ATI  Vf:      Iffji 

SUM    =    su"    f    ,--!;(-•;)  *J.<(:<> 
'       Gu'^iA"^r":;i!"->.r  V"  ^ukT 

F,UVAe     =     SiJ.-'AO     +     A:i3i   -•'(<)  ) 
Su:PCr    =     SO'-J'tT     t     PCI--->(<) 
~C *  i  it" i^rlv LT:TiL  aT' ■    ~il-r6 ~"S"o I ■  ~ ~Ji^"' ~T"h E" ~ 0~ J "i)zuvX f  ;  r:' N  J 
suf-ix  =  suy  <t  a:  {  <  ) 

C    U.^DAF:     the    CO\ST\:iTS 


U044--i  [  =  1  .  \ 

4  4  'J    C  I  E  (  I  )  =  0  .  y 
i:ij  i'i4  !  =  I  ,   , 

9  4  4     C  I  L  (   n  =CI  L  (  I  )  t  T  T  ■■  (  I  .  J  )     C  !  t 
DO'32;.  !  =  1  .'I 


EXPOSm  (cont'd. 


CI  (  1  )     =    CIU(  I  ) 

CONTI-NlJt-: 

E'^y     -     SUMA/ND 


H(  :  ,  1  ).!■.£■!(  K) 
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SU.".A.'3CTeH(K)  ) 


ASCKH  =  SUVAil/fU) 
PEk'M  =  SU.MPCr/Ni) 
00    20         1  =  1.   40 


VAKE   =   s=r)So/r;ii 

SOEVE     =     S()C!  r(  VAJ~  1 
"SQ.<trsUM'/ND) 
!) 


SDgV 

V)i-II  TE  (  J,  lOo  )  (  SUV) 
'  C^*  *V* C'^L.OTirS'TT~f  Ho" "  V  Vi T  A N C e 
VA'i=SUV/  KD 
SDFST     =      SOJTlVAf^) 


QF     THS    FUi^ECiVSTS 


»tMTF(    i.  iiP.  )  VAi( 
WRITE  (.1.1)         301-- ST 
WH  I  TE ( J,? )        EHM 

'AIM  TE(   3.  1) 
W^ITEIS.A) 
ViRITEC  J.'3) 

AtlEHN' 

VAi<E 

SDEVE 

■aHITE(   J.O  ) 
WHITCI 3.7) 

■^ETUK'.M 

SDEV 

END 

E3CP0SM0   (opnt'd.)„ 


5 _^ 200_ 

SUoi"U;UT  !  -Ih    PUoi  (NO,  A,N.:'  ,nL  ■f-.i:  j 

l)I"CMnlO  -J     OU  r  (  1  on  .  Yl'  ■)(  1  1  )  ,  ANG(    U  .  AC  n 

c 

;>     >^0!".M(1H      .  F  1  1  .  'I  ,  i.il+  .lOlAl) 

"!   r-c  r-i,r  (  m-  ) 


5    KOJMATC lOA I ) 

7    FLT-'MAT(ti(     ,l6:<.101h  + + + 1- 


l + + + : ^ + <■  ) 

H    FmHAT(  IMC  ,  '1  <  ,  I  ir-  1  0.  4  ) 

'I    FU^.MaT     (1H     ,  [hit  ,  1  o.<  .  1  0  I  a  1  ) 


c 

c  

c 
c 

I F  (  '  1 5  )      16,      I  'i  ,      10 

c 

C  SGK'T     BA"")"     v/Ar-lI  AfcLI     OAIA     IN     ASC^IMUKG     Ci-iOfER 

C 


10   oc;    10    1  =  1  .^; 

DCJ     l'^     J=I  ,'■ 

iF(A(n-\(j))    I'*.    I'l.ii 
TT~U=  I-  1 
LL=J-^; 
DO    12    K=i,,; 


L  =  l_  +  -J 
L.L-LL^^; 
r-  =  A  (I. ) 


A  ( 1.  )  n''  t  LL  ) 

12     A(LL)=F 

I'l     COJTl'lUi^ 
^._„  _^._.-._._.^ 

C 

C  Ti_;ST     NLL 

~C 

lo     iri-lLL)     P.O.      l>i.£0 

c  i-'.;rM    riTLc 

c _ __ 

^-,— --J-J  fT?-" ("57  \~)Tvi  ' 

c 

C  OKVi-LCP     dLaD:;     AN^)     digit:     r-Z.i     P'3I\riK3 

r-%£.'llNO       l.' 

utit  T'     (  111  ,  J  ) 

i^CAU     (  1  0  ■  i  )         t)LA^:K  ,  (  AfJi;  (  !  )  ,  I  -  I.  .  ^) ) 


c 

C  Fi-;r.    re ','_':   fc-!    ja'sc   va;  i  •■jl-: 

' ~~  /.'ryz'^C-T'iO'n-^Ti  i" l/fr T:  \r(  ji.l-~i  ) T" 

c 

c  Finn    GCAL-    f-t:"   c»(, s ■'.-.; a -■  i '.eLK; 
"       '    EXPOSMO    (ocnt'd.) 
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V-  AX  =  Y,VI  :j 

'.Tl     '40     J=  ■'!  ,  y2 

IF( '\(  J  )-Y."  [;-.)    ■da,i'>,26 


ii:   in .\( .1  )-y;- Ail   "»i'o7<roT3"d"~ 

L'(      YMIN  =  A(  J) 

GO     TO     AO 
^0~YT.'.\X'HaTjT~: 
'lO    CU-!ri:JUE 

Y5C^L=(  Y'l'iX-Y.vi.'i)/!  CO.O 


C  FIND    fiASti     VA«IA'3LS     P  i?  I  ^1 T     PCSITICN 

c 


XF_<  =  A(   Jl" 
(1Y  =  i1-  1 


1  =  1 

i  a   r  =  I  - 1 

XiJ^'^XH  +  FtXSCAL 


~rr  rA  (  L  J  -X.MT"  'FO  ,  i.  -5  ,  TCT 

c 

C  FIND    C,'GES-VAV  I  AciLES 


A<)   l.  s\.  =  o 

so     IF     (L-N)     :3I.5  1,tiO 


'ji   DO  '--i    i/=i,uri 

Sfj     IHJK    1  >;)  =:.LA^JK 
■37    DO     60     J=  t  ,  ^Y 

JP=(  (  A(LL)-Yt.   IN  j/YOCADt  I  .  0 
GUr(  JI')  =  A'JG(  J) 

L=L+I 


c 

ir(LSV.)     62.Lil,0)2 

"'l^i"^TT^rr.iVi:l~XpTr,TciTTrf7T7lZ^r,TTnT" 
Lr.w=i 

rF(  a(l)->:?-!-xS(;al  )    30,^o,tiO 

IF  (  Ad.  J-XPJ-xaCAL}     'jO.  "0,dC 


7U"  "^"rnT-Tr3V3") 

fiO     1  =  1  +  1 

IK(  I-^JLL)    .-ts,.;/*  .a-i 
"Tv*^  X .'  ■)  =  M  IT) 
GO     TO     -iT 


c 

pt' 

l-,T     C(C 

c 

tic    -  r<  1 

[  rhi 

:  1.7) 

' 

YP. 

■(  1  : 

l  =  Yy.  IN 

on 

90 

<M=1 .Q 

'  c  ?.  ■;  -  V  A  >i  f  A  H  L !  ■  S     J  U  '■'  L-;  -I  - 


90     YP'(KNM)=YP->(i^\)4-YSCAL*I0.0 
EXPOSVO    (cont'd.) 
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YPi^l  1  1  )=y\:\x  _  

■■■tiiru-!--! 


_EX.PCS;iP_  (opnt  •  d_il . 


DI..'|":ns  ILf.       ;(.<  (  700  )  .  VY  (  1   it' )  ,  xr- ^<   :■   ),YrY(]) 203 

cu'rr,r.ic\   xf  Y(o  ) ,  ji- ',  M  c  ) ,  j.w:;  ( •. )  ,  mV.  )  ,!.'i;x^:V  (  f ) 

DIf-EKCU'N     Yh.i  r  i  ir.C  )  .PLT  (  V;G  ) 

_  _DU--E\s\c-.   ^(icc),  i(  1 1  .■^ )  .  '(  1  cr  )_    

cir'e'vE  i(;\    itz:,i::\  i  vc  )  .px(  ioc  )  .i-:^/.  i  tci.  )  ,r-xxx  (  2(, Cj  ) 

I    r  c  'w,,  r  (  1  6. 1  J  ) 

i'   ri^^;.'M     (;-b.O) ■ J 

J   rcji'i- /\  r !  iH'vc'. '<,:'('!-',■>:.  xLv; :  ■".    cj-    v.\  u  ;-.nc;.' ) 

«  r  cj-',';».t  (  iHC  ,  3  <,  i<;hiiji.  Kii  j'^  v  . -j '.  it  it.i.  .7x,;Hr:Er,u!;-  uf  .  7v  .  t-t  ru;-  (k. 
_  1  .''H.'-'rtAN .  lox  .  5^i2£/  jjx  ,  7;iF;;'7.:nL;- .'^x .  ''hscu  -I'.i-S .  /•  < .  enrcui.'"'^  ■  vx .  Shv..i.t 

5    Fpflf  ATI  J'CH.I.'        nLC     TIJ     UFG'   ; /J":  I  I  ,f.  .  1  i^<  .  !  i:  .  F  1  '  .  S  i  F  1  7  .  J  .  F  1  -  •  S  ) 
e    rU3hM(3?H  D=VI/iTIi.\     ^.   CUT     M'-T.-ir':'^.  IC\  ,  I6.F  1  7.b  ,f"  I  7.0/// ) 

7    fC-i.'-'nTC  inC  ,  4h-'J(:     =  ,K  1  0  . 'j  .  '  CX  .tl-^^i^.F  IT  .5  .  l_t'Xl' '•  Ij^L"  T     -  ■  F  1  (    .  '.'  .  1  C  ■<  ■  '■  > 
,  5  ,  1  C  X  , 'IH!,',     =,FlC.-;/Cl-     i;  ;     =,flO.=;j 

(  i!-.i-/.j-4x  ,;' ;ii.,,£cc  ;.=i;k  I  T  It:  V    of-    -  =  situAL/) 

(jho.  jx,  I'.hsci;  <cr    ^r-^  ^/■^  :  [  \  n  l:..  .  7  < .  s.i-n':'3  iEf    ^^  .7x  ,':nsi.'   ur 

1  X  .ih'VCAN.  1  ex  ,  1HF/3.1X  .  /f'KIJcHF.rjC ..  .  -iX  ,  ?HSCUai.  =  S  t  7X  .  tl-SCUAVfc  ,  /x.  thVAL 
2) 
le_Fr;  ■<.■■' AT     (I4H0_      lack     CF     FIT.lf--;  ■lc.FI^.5.F17.'i.F13.b) 

".^•cp.  i-'x  ,  i  t  ,r  I '.b  .F  1  7.5//) 

GF.TlN'.Ti:  "FS  1  CU'-L/ ) 


1 

1     -, 

,  F  1  C  . 

le 

1  (jT.-  M 

17 

ft- 

'r  AT 

I<:j 

f  C^l'  '■  T  (  1  3H            PbHC     ,";i 

if^'l 

?C 

FC:«r.'AT    (  7F  I  ^.s) 

il 

F  U^'f'AT     (  T<.l'               ACTL. 

■iL 

2r 

FC:'-i!-AT      (F?r,./-.) 

ct 

■VTAO      (1.1)      N  ,KC 

1  F  (  ^.  )     ?  Cj4  .  ?■  C  •i  .  9  5_ 
l^F'i^D     (  1  ,J)      (  VY(  I  )  i 

nr,   A'.-n    :  =  !  .k 
xx( I  )    =    1 .  c 


XX( I+N)     -     I 

XI      =      I 

XX(  I  +2*t.<)      -     X  !<■*? 

XXC I+3*N)      =     XI**3 

XX(It4*K)      =      X£**A 

aoc 

XX(   I+SSN)      =     X I**b 

CALL     GTP.iD     <  XX.  XX  ,  XPX  .\  ,o  .C  ) 

CALL     GTP^^Il     C  YY.  YY  ,YPY  ,N,  I  ,  1  ) 

CALL     GTI^-VD     (  XX  .  YY.yPY  ,N.6  ,  I  ) 

CALL      '■lll'j      (XPX.fj.DCT      ,  JV.  ■;  1  .  J   ,V;    ■ 

CALL     GTP,'0     (  XPX  ,  XI'Y  .r  .:  ,6  .  1  ) 

CALL     GTP.'O     (  ■_.  .XPY.QPXPY.o,  1  .  1  ) 

YSIJS;    =     O.C 

DC,     IOC     K=l  ,N 

ICO 

Y  S  U  !■■:     =     Y  =  U  I'     t     Y  V  (  K  ) 

IN     =     \ 

CFAC      =      ( YSUy**2 ) /TN 
_Cr.SDK_=     'JPXPY(J|       -CFAC 
.^G>-S     -     CS'iO'J 
CS'-.T!.     =     YPY  (  1  )        -JPXPY(l) 

ijL'';^s,    =    c-;r^iL/(^TN-6) 

J-''l     -      5 

J'VL     =     N-6 

F  V  IL     =     ^  : !'  S  /  ;£U^  S 

V  i'<  I  r ::    (3.3) 


Wi-'IT'^      1 

:3.<*i 

>>/iT':    1 

;  3 ,  =i ; 

1      J  -1 G  .  C  J  S  0  -! 

n  ►' !  T  :    1 

:  3.  J : 

!      .)-.'L  .CiG -(L 

^  UT  :    ( 

:  i  ■  ."■ ; 

1      r  (  1  )          .(,( 

ai'j   ?co 

1  =  1, 

.  \ 

(   *  )  ,L  (  j)  .;:(  t>) 


pb3gru'  no.  f-5 
i;atr2Gs 
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It  ■^^  I T*'.'  (  3  M  X .;  (  5  *  N  +  [  )  *  iJ  ( ■'! ) 

uii    301     J  =  I  .w 
J'Ol_-'LT(J)      =      'M.~±J) 

DC     2C?     --^l.K 
3C2    Pl-T  (  ><  *  V  )     =     Y  Y  (  i<  ) 

mi     ?C3    L=l  .\       ^ 

~2C3    .^LT(L  *;>*■.  )      =     VI-ATCl) 

c  XL,-    'l;  r    (  N :  ."LT  ,\  ,  )  ,:>  3, 1  ) 

or    -TCO    1  =  !.N  


3t;C     '■!;SIO(  I  )iYY(  I  )      -     YHATCn 

V  «  Tin  !  ? . ;  2  r   (  .'Ts  I  C  (  1)  ■     I  =  1  .  \  ) 
GC    TC    ';s 
_acj4_  sjnp     


c 

C  SUP.^CUT  I  Ki;    GTP«0 


c 

C  PL.'^rr'SH 

.C  PPr->LLT  IPLY     4     jirKEi^AL  JwlJQ^iJ<_  _ii'^  Jj 


C  GtiVE.;  -iL     .'    AT'tl  J< 

c 

_C UJ^AiiE ^ 

C  CALL     GTP>iC(  'l.a  .!■' .Is  •!.' .L) 

c 

C  DhSC^IPTIl.K     CF     P'K,tl\r  =  TZfiS 

C  a"  -~  KA.vr^CF   >  I  kET      INPUT   ".'-'ATPIX 

C  E.     -     KAIE     CF     SECCM/     IKPLT     ^^TPIX 

C K     -     NAr  h_    CF     (JLTPLT    J- M-!_I  X 

C  N     -     N  L  !■■  h:  c  R     OF     P  P  V-  S      IN      l^     »U.  C     e 

C  h'    -     KUreEF     (JF     CCLLINS     IN     A     AND     SCv 

c  L  -  ml;  f.'Ci;  LF  ccll'-k:  in  ~    anc  r 


c 

■  C  R  e  V  iMi  K  s 

c  wATPix  R  cannct  he  in  if,"  sa;'=  lccaticn  as  NAr-^lX A_ 


C  N'AT^atX     R     CANNCT     f£     IN     The     SAf-L    LCCATIC\     AS     ;.lAI-iIX     F 

C  ALL     VATKICES     ycST     Be     SICPci:     AS     C  =  NFftAL     r-^ATfilCES 

c  


C  SLbPOLTI.N^S     AKD     FLNCTIl.N     SLbP -^i;.' ^  AI' E     =;fCLISlI; 

C  NLNc 

c  


C  KETHUO 

C  ^•ATPIX     T^ANSPtr.';     or      A     IS     NLT     ACTUALLY     CALCULATFD.      IN'Fr';AC 

■C  FLFf/^NTS     CF     yAT-^IX     A     A.^f     TAKirv     CClIVNaISF     kATFF.? THAN 

C  ^icv7.\  I  EC     F  C"~PcVf7/ULT  i  PL  ICAl  !CN     [- Y     :   A  T  J  I  X     E. 

c 

C      ^ 

SLe'JrLl  I  N-Z     GTPi^  C  (  A  ,  t  ,  ^  ,M  .  .'•■  .L  ) 
Lll'--E-<Sn  \     AJ  I  CO  )  ,J-  (  1  fiO)  ,  V  (  1  cc  ) 

MATRSGH  (cont'd.) 
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1^=0 
_       Il=-N . 

1  J=C 

' ^ '  •'^z.'.'■i^■'■^ 

DCi     1  c;     .J=l  ,1' 

■ Ji^L^-*J _ 

H  {  I  p  )  =  c         ' 

nu    1 c    [=1 ,N 
ij=ij+i 

IC    i^(  !H  )  =  5(  I  n +«(  t  J)  tii(  la  ) 

^(GnjKj^    

FND  . 

c    ' .". . ,.;, 

c 

c ?_'^'i ''f'-dL^-^i-  "jhy_^ 

C  '   '       -         - 

C  PLWOSP 

C  INVI4TAyATtiI.< 


c 

C  USAT£ 

_C C ALL^   t'  IJ\ V_(  A  tj^tp  iL  ,^J 

C  '  "^  '""  '"  ' 

C  DESC>^IPTIC^    CF    P^■<A:^t.rz■■iS 

.c ^--"--'-'^■^y--'^''*. !:'"'_•., ,p'ij_f  <cvci:    in    ccyp'jTAriCN   anc   ^c.pl,\c-c   i 

C  «£SLLTAi\T     IKVE^it.                   "                                           '                      ^ 

e  N     -     C !-!  0  e :  i     C  F     .■•  A  T  R  I  <     A 

_C i)_  -  J? e  S LL  T  AM     Oi  X li  ) .''_ I  N /IN  1^ 

C  1.    -     HL^K     J-Crci     CF     L'f.N.GTH     N                               1 

C  y     -     v.CRK     viCrCR     CF    LtNGTh     \ 

_C 

C  RS.VARKS                                                                                                                                                                          ~ 

C  VAT^ilX      A     yoST     f?     A     Gi=\';-^A|.     VUT)-fI>! 

c 


C  NC\;£ 

jC ; ■_ 

C  M'iTMCD  " 

C  THf!     STANCI-JU     GAUSS- JC«L'A\     .vzrhCC     IS     UTfiC.     THE     L't  T'r-  .' I  \  AM 

_£ J  3     ALSC   j:ALi^L\TeO^    A_  i;;:rC^INJVNr     CF     iF-'C      INCICATSS     Tt-AT 

i:  TH£     yMiilx     13     SIK.luLA'J^ 

c 

c  , 


c 


D_If'E\S  ir.N     A  (  1  Co  )  ■!_(  ^  )  ..•'  (2  ) 


:C  IF     A     Dfn.dLtr     PDECISICN     VFJSl'.N     C  f"  "r  F  f  3 "  -^C  U  f  I  NT     Is"cFrrrPSc\     TF 

C  C     PJ     COH...N     1     SHCLLO     iiC      -Ili.v.-y^D     F,;L.y     THE     CCUFLE     P-'ECISICN 

.C  __  STATEVFNT     v.HlCH    FCLLQaS. 

KATRESR  (cont'd.) 
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c 

c     ccunLE  p-;ecisu:n  a,d.hig\  ,h:lc 

jC 

C  Th£  ■:  rusr  AL'n  h;  Pr-frM\-i:     ri-CK  CCLrLL  P-hCISin\  ET  ATrl-'C.  rj 

C  APPEARING  IN  CTUti'^  "CU^1^^S  tSCC   I  ,N  CCNJLNCTiCI^  V  I  Th  ThIS 

C  RCUTIKE. 


c 

C  THE     DHUBLJ    PKfcCIPlCK     VLHEICN     CF     THIS     :>UhSCL,TIkE     f-L'ST     ALPC 

_C CCNTAIN     DCXbLC     PREJCIPIC^K     FO-;  T  ^J  AN     FL.\C  TJj;  NS  .        AbS      IN     STATEvEN' 

C  10     VUSr     BE     CI^ANGEC     TO     UAHS. 

c 

c 


c 

C  SFAPCh     FC'^     LAfJoeST    rLEf-'ENT 

c 


0=1  .tl 

NK=-N 

^ DO    80    K=1.N 

_NK=NKt\ 

L(kV=k 

^,(K)  =K 

KK=Nt;tK 

e'lGA  =  A"{K<) 
DO    2C     J=<,N 

IZ  =  N*(J-I) ._ 

DC  2C  I=K,N 
IJ=I Zf [ 
■   IC  IF(   AOSdUGA)-  ACSCACIjn)   i; 
'  13  13  IGA=4"(  IJ) 
L(K)=I 

^^(_K  )=J 

20  CONTINUE 

c 

C JJ\TEJ<ChANGE     CCftS 

C  '        •    ' 

J=L(K) 

IF ( J-K )     3b . 35.25 

iS     KIsK-N 

DO     30     1  =  1  ,N 

■ K  I=KI  tN 

VcLD=-A(Vn 
Jl=Kl-K<J 

A  ( K I )  =  A  ( J  n 

30     A  ("JI  )    VhCLD 

c 

:  C IK  T  li  1^  Oh  a  n  G  E_  C  C  LL  ;•■  .N  S 

c      ' 

35    l=;^(<) 

IF ( I-< )     H5, 4 i , 32 

3a     JP  =  Nt (  I- 1  ) 
DC     40     J=  1  ,K 
JK^NKtJ 
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Jl=jp+j 

HCLD--A(JK) 

.1  (  JK  )=A(  JI  ) _ 

10      I ( JI )      =hCLb 

c 

C DlVlCf     CCLL"K     SY     "'I^IJS     PIVCT     (VALU?     CF     FIVCT     ELE^EnT      IS 

C  CCNTAINCD  IN "f IGA ) 

■c 

_5_5     lF(eiG_A|    tb.'ib,'iS_      __ 

46    D=0.0 

Jte     DC     55     I_=1,N ' 

IF(  I-K  j     50,33  .'oO^  "  '  [ 

5C     IK=NKtI 
A(  IK  )JA  (  IK  )  /  (  -C  IGA  ) 

SS    CCNTINUe  "■ 

c 

C n  E  0  UCE     F  AT  R  IX 

C     "  -  - 


DO     65     1=1 ,N 
IK=NK«I  _ 

H'GLb=A(  IlO 
I.J  =  I-N 
00     65     J=l  ,N 
IJ=I J+N 


lr(!-K)     6  0,65,60 
CL      IF( J-K )     62,6  5,62 

i-JJi'J^  I  J- J  +K  _ 

A(  I  J  )=HCLO*A{  KJ  r+M  I  J) 
fS     CQNTINLE 
C^ _ 

c  D I V  i  DE   '^ c  A   ij  Y   pTVc  T 

c 

J<J  =  I1-N 

00     75     J=l,"i.i 
K J=K J+N 

n^(  J-<^      70.  7S,yO 

7r  A(K.J  )  =A  (  K  J)/bICA 
75  CCNTIN-LE 
C 


c        PHCcucr  cr  pivcts 

c 

_D=^DJ'f_I_GA^  _    

c  " "  '  '    '"" 

C  REPLACE  PIVCT  UY  ;!ec  I  P.TC  CiL 

C 

A  (  K  K-  )  =  1  .  0  /  h'  I  G  A 
CO  CnMTIMUh 
c 

c  FINAL   ■<cn    AND   c c L u y frTrrTtfi^cTi A XgI; 

c 

K=N 


MATREGR  (cont'd.) 


208 


IOC     1-  =(K-1  ) 

1  F  (  K  )      1  'j  G  ■  1  3  I.  .  1  t 
J  Ci     I=L( f  ) 

I  F  (  I  -  ".  )      1  2  0  .  1  Z  t  . 

jr7  =  N*t  l-l  ) 

ro    lie    j=i fN 

JK= JC+ J 

hUI.O  =  *  (  JK  ) 

J  I  =JK  +  J 

A( JK)=-A(  J  I  ) 

1 1 0 _ A  { J 1 1  =1  inLc 

120     J=v( K  ) 

1F( J-K  )      1  CO ,  1 CC , 
liibK  I=K--N 

DO     130     I =1 .N 

K  I  =  K  I  +  N 
IIGLO  =  A(K[  )       


A(<I  )=-A(  JI  ) 

nr.     A(  J  1  )  -HOLD 

GC  TO  100 

150  RETUr^N 

END 

MATREGR  (cont'd.) 
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The  ability  to  forecast  with  reasonable  accuracy  the  future  demand 
for  an  item  of  production  has  value  anywhere  that  an  inventory  must  be 
maintained.   This  research  investigates  one  approach  to  the  problem  of 
forecasting,  general  exponential  smoothing. 

The  logical  basis  for  forecasting  the  future  demand  of  an  item  is 
the  demand  behavior  of  this  item  in  the  past.   Such  past  behavior  may 
be  summarized  in  the  form  of  a  time  series,  a  sequence  of  observations 
of  demand  taken  at  specific  and  equal  intervals  of  time.   We  assume 
that  this  time  series  is  composed  of  a  systematic  component,  which  may 
be  described  by  a  linear  combination  of  deterministic  functions  of  time, 
and  random  variation.   The  deterministic  functions  of  time  are  known  as 
fitting  functions,  and  their  sum,  along  with  the  appropriate  fitting 
function  coefficients,  makes  up  the  model  of  the  time  series. 

Assuming  that  recent  data  is  more  representative  of  the  future 
time  series,  we  wish  to  weight  the  data,  giving  the  most  recent  obser- 
vations more  weight  than  older  data.   We  do  this  by  assigning  a  weight 

It 
a6  to  an  observation  made  k  sampling  intervals  before  the  present  time, 

where  ex  and  6  are  constants  whose  sum  is  unity.   Now,  we  may  control 
the  general  exponential  smoothing  response  rate  by  merely  varying  B, 
the  discount  factor.   When  6  is  smaller,  we  give  recent  data  more 
weight,  and  the  forecasting  system  may  respond  more  quickly  to  changes 
in  the  actual  data. 

General  exponential  smoothing  fits  a  least  squares  curve  to  the 
time  series  much  like  statistical  regression.   It  does  this  by  determin- 
ing the  fitting  function  coefficients  in  the  model  that  will  minimize 


the  sum  of  weighted  squared  residuals.   The  main  differences  between 
these  two  methods  are  that  general  exponential  smoothing  operates  on 
weighted  data,  and  the  fitting  function  coefficients  are  re-computed  each 
time  new  data  becomes  available.   Hence,  this  process  is  able  to  adapt 
to  small  changes  in  the  true  time  series  process. 

To  simplify  the  computations  involved  in  making  a  forecast,  the 
fitting  functions  must  consist  of:   (1)  trigonometric  functions,  (2)  poly- 
nomial functions,  (3)  some  exponential  functions,  (4)  mathematical  combin- 
ations of  the  first  three  types.   In  addition  to  this  restriction,  the 
number  of  time  series  observations  available  must  be  at  least  equal  to 
the  number  of  fitting  function  coefficients  to  be  computed. 

General  exponential  smoothing  offers  several  advantages  over  other 
fore-casting  methods:   (1)  Response  rate  is  easily  controlled;  (2)  All 
past  data  is  considered  and  contained  in  one  word  of  information;  (3) 
Model  coefficients  are  re-computed  at  each  sampling  interval. 

The  data  used  herein  consists  of  four  time  series  of  demand  for 
plants  produced  and  sold  by  Manhattan  Wholesale  and  Retail  Floral  Companies. 
For  each  time  series  various  combinations  of  models  and  discount  factors 
were  tried  in  an  attempt  to  determine  the  combination  that  best  forecasts 
the  time  series.   Forecasts  were  evaluated  on  the  basis  of  three  efficiency 
criteria:   standard  deviation  of  the  forecasts,  mean  absolute  error,  and 
standard  deviation  of  the  error.   The  effects  of  model  and  discount  fac- 
tor choice  on  fore-casting  efficiency  were  also  studied  in  an  attempt  to 
gain  a  clearer  understanding  of  the  general  exponential  smoothing  techni- 
que. 


^ 


All  computations  are  performed  on  the  IBM  System  350/50  digital 
computer  at  Kansas  State  University. 


